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1 Introduction 

Many properties of condensed matter systems can be 
calculated from solutions of the stationary Schrodinger 
equation describing interacting ions and electrons. The 
grand challenge of solving the Schrodinger equation 
has been around from the dawn of quantum mechanics 
and remains at the forefront of the condensed matter 
physics today and, undoubtedly, for many decades to 
come. 



The task of solving the Schrodinger equation for 
systems of electrons and ions, and predicting the quan- 
tities of interest such as cohesion and binding energies, 
electronic gaps, crystal structures, variety of magnetic 
phases or formation of quantum condensates is noth- 
ing short of formidable. Paul Dirac recognized this 
state of affairs already in 1929: "The general theory of 
quantum mechanics is now almost complete . . . The 
underlying physical laws necessary for the mathemat- 
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ical theory of a large part of physics and the whole 
chemistry are thus completely known, and the difficulty 
is only that the exact application of these laws leads 
to equations much too complicated to be soluble." [1] 
Arguably, this is the most fundamental approach to 
the physics of condensed matter: Applications of the 
rigorous quantum laws to models that are as close to 
reality as currently feasible. 

The goal of finding accurate solutions for stationary 
quantum states is hampered by a number of difficulties 
inherent to many-body quantum systems: 

(i) Even moderately sized model systems contain 
anywhere between tens to thousands of quantum 
particles. Moreover, we are often interested in 
expectation values in the thermodynamic limit 
that is usually reached by extrapolations from 
finite sizes. Such procedures typically require 
detailed information about the scaling of the 
quantities of interest with the system size. 

(ii) Quantum particles interact and the interactions 
affect the nature of quantum states. In many 
cases, the influence is profound. 

(iii) The solutions have to conform to quantum sym- 
metries such as the fermionic antisymmetry 
linked to the Pauli exclusion principle. This is 
a fundamental departure from classical systems 
and poses different challenges which call for new 
analytical ideas and computational strategies. 

(iv) For meaningful comparisons with experiments, 
the required accuracy is exceedingly high, espe- 
cially when comparing with precise data from 
spectroscopic and low-temperature studies. 

In the past, the most successful approaches to ad- 
dress these challenges were based mostly on reduction- 
ist ideas. The problem is divided into the dominant 
effects, which are treated explicitly, and the rest, which 
is then dealt with by approximate methods based on 
variety of analytical tools: perturbation expansions, 
mean-field methods, approximate transformations to 
known solutions, and so on. The reductionist ap- 
proaches have been gradually developed into a high 
level of sophistication and despite their limitations, 
they are still the most commonly used strategies in 
many-body physics. 

The progress in computer technology has opened a 
new avenue for studies of quantum (and many other) 
problems and has enabled researchers to obtain re- 
sults beyond the scope of analytic many-body theories. 
The performance of current large computers makes 
computational investigations of many-body quantum 



systems viable, allowing predictions that would be dif- 
ficult or impossible to make otherwise. The quantum 
Monte Carlo (QMC) methods described in this review 
provide an interesting illustration of what is currently 
possible and how much the computational methods 
can enrich and make more precise our understanding 
of the quantum world. 

Some of the ideas used in QMC methods go back 
to the times before the invention of electronic comput- 
ers. Already in 1930s Enrico Fermi noticed similarities 
between the imaginary time Schrodinger equation and 
the laws governing stochastic processes in statistical 
mechanics. In addition, based on memories of his col- 
laborator Emilio Segre, Fermi also envisioned stochas- 
tic methodologies for solving the Schrodinger equa- 
tion, which were very similar to concepts developed 
decades later. These Fermi's ideas were acknowledged 
by Metropolis and Ulam in a paper from 1949 [2], 
where they outlined a stochastic approach to solv- 
ing various physical problems and discussed merits of 
"modern" computers for its implementation. In fact, 
this group of scientists at the Los Alamos National Lab- 
oratory attempted to calculate the hydrogen molecule 
by a simple version of QMC in the early 1950s, around 
the same time when a pioneering work on the first 
Monte Carlo study of classical systems was published 
by Metropolis and coworkers [3]. In the late 1950s, 
Kalos initiated development of QMC simulations and 
methodologies for few-particle systems and laid down 
the statistical and mathematical foundations of the 
Green's function Monte Carlo method [TJ. Eventually, 
simulations of large many-particle systems became 
practicable as well. First came studies of bosonic 
fluids modelling 4 He jSHT], and later followed investi- 
gations of extended fermionic systems exemplified by 
liquid 3 He [SI [S] and by the homogeneous electron gas 
[TU1 ITT] . Besides these applications to condensed mat- 
ter, essentially the same methods were in mid-seventies 
introduced in quantum chemistry to study small molec- 
ular systems [12] [13]. To date, various QMC methods 
were developed and applied to the electronic structure 
of atoms, molecules and solids, to quantum lattice 
models, as well as to nuclear and other systems with 
contributions from many scientists. 

The term "quantum Monte Carlo" covers several 
related stochastic methodologies adapted to determine 
ground-state, excited-state or finite-temperature equi- 
librium properties of a variety of quantum systems. 
The word "quantum" is important since QMC ap- 
proaches differ significantly from Monte Carlo methods 
for classical systems. For an overview of the latter see 
for instance [14]. QMC is not only a computational 
tool for large-scale problems, but it also encompasses a 
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substantial amount of analytical work needed to make 
such calculations feasible. QMC simulations often uti- 
lize results of the more traditional electronic structure 
methods in order to increase efficiency of the calcu- 
lations. These ingredients are combined to optimally 
balance the computational cost with achieved accuracy. 
The key point for gaining new insights is an appro- 
priate analysis of the quantum states and associated 
many-body effects. It is typically approached itera- 
tively: Simulations indicate the gaps in understanding 
of the physics, closing these gaps is subsequently at- 
tempted and the improvements are assessed in the next 
round. Such a process involves construction of zero- 
or first-order approximations for the desired quantum 
states, incorporation of new analytical insights, and 
development of new numerical algorithms. 

QMC methods inherently incorporate several types 
of internal checks, and many of the algorithms used 
possess various rigorous bounds, such as the varia- 
tional property of the total energy. Nevertheless, the 
coding and numerical aspects of the simulations are 
not entirely error-proof and the obtained results should 
be verified independently. Indeed, it is a part of the 
modern computational-science practice that several 
groups revisit the same problem with independent 
software packages and confirm or challenge the results. 
"Biodiversity" of the available QMC codes on the scien- 
tific market (including QWalk [15], QMCPACK [32], 
CHAMP p7], CASINO pi], QMcBeaver [19] and oth- 
ers) provides the important alternatives to verify the 
algorithms and their implementations. This is clearly a 
rather labourious, slow and tedious process, neverthe- 
less, experience shows that independently calculated 
results and predictions eventually reach a consensus 
and such verified data become widely used standards. 

In this overview we present QMC methods that 
solve the stationary Schrodinger equation for con- 
densed systems of interacting fermions in continuous 
space. Conceptually very straightforward is the varia- 
tional Monte Carlo (VMC) method, which builds on 
explicit construction of trial (variational) wave func- 
tions using stochastic integration and parameter opti- 
mization techniques. More advanced approaches rep- 
resented by the diffusion Monte Carlo (DMC) method 
are based on projection operators that find the ground 
state within a given symmetry class. Practical versions 
of the DMC method for a large number of particles 
require dealing with the well-known fermion sign prob- 
lem originating in the antisymmetry of the fermionic 
wave functions. The most commonly used approach to 
overcome this fundamental obstacle is the fixed-node 
approximation. This approximation introduces the 
so-called fixed-node error, which appears to be the key 



limiting factor in further increase in accuracy. As we 
will see in section [5] the fixed-node error is typically 
rather small and does not hinder calculation of robust 
quantities such as cohesion, electronic gaps, optical ex- 
citations, defect energies or potential barriers between 
structural conformations. By robust we mean quanti- 
ties which are of the order of tenths of an electronvolt 
to several electronvolts. Nevertheless, the fixed-node 
errors can bias results for more subtle phenomena, 
such as magnetic ordering or effects related to super- 
conductivity. Development of strategies to alleviate 
such biases is an active area of research. 

Fixed-node DMC simulations are computationally 
rather demanding when compared to the mainstream 
electronic structure methods that rely on mean-field 
treatment of electron-electron interactions. On the 
other hand, QMC calculations can provide unique in- 
sights into the nature of quantum phenomena and 
can verify many theoretical ideas. As such, they can 
produce not only accurate numbers but also new un- 
derstanding. Indeed, QMC methodology is very much 
an example of "it from bit" paradigm, alongside, for 
example, the substantial computational efforts in quan- 
tum chromodynamics, which not only predict hadron 
masses but also contribute to the validation of the 
fundamental theory [2Ql [21] . Just a few decades ago it 
was difficult to imagine that one would be able to solve 
the Schrodinger equation for hundreds of electrons by 
means of an explicit construction of the many-body 
wave function. Today, such calculations are feasible 
using available computational resources. At the same 
time, there remains more to be done to make the 
methods more insightful, more efficient, and their ap- 
plication less labourious. We hope that this review 
will contribute to the growing interest in this rapidly 
developing field of research. 

The review is organized as follows: The remainder 
of this section provides mostly definitions and nota- 
tions. Section [2] follows with description of the VMC 
and DMC methods. The strategies for calculation of 
quantities in the thermodynamic limit are presented 
in section [3] Section [4] introduces currently used forms 
of the trial wave functions and their recently devel- 
oped generalizations. The overview of applications 
presented in section [5] is focused on QMC calculations 
of a variety of solids and related topics. 

1.1 Many-body stationary Schrodinger 
equation 

Let us consider a system of quantum particles such as 
electrons and ions interacting via Coulomb potentials. 
Since the masses of nuclei and electrons differ by three 
orders of magnitude or more, the problem can be sim- 



3 



plified with the aid of the Born-Oppenheimer approxi- 
mation, which separates the electronic degrees of free- 
dom from the slowly moving system of ions. The elec- 
tronic part of the non-relativistic Born-Oppenheimer 
hamiltonian is given by 



Zj 



Ti - Xi\ 



E 



i 



(i) 



where i and j label the electrons and I runs over the 
ions with charges Zi. Throughout the review we em- 
ploy the atomic units, m c = h = e = Atteq = 1, where 
m c is the electron mass, — e is the electron charge 
and eo is the permittivity of a vacuum. We are inter- 
ested in eigenstates of the stationary Schrodinger 
equation 

H\V n ) =E n \V n ) . (2) 

Colloquially, we call such solutions (either exact or 
approximate) and derived properties collectively the 
electronic structure. 

An important step forward in calculations of the 
eigenstates was made by Hartree [55] and Fock [53] by 
establishing the simplest antisymmetric wave functions 
and by formulating the Hartree-Fock (HF) theory, 
which correctly takes into account the Pauli exclu- 
sion principle [5~C 25.. The HF theory replaces the 
hard problem of many interacting electrons with a 
system of independent particles in an effective, self- 
consistent field. The theory was further developed by 
Slater [55] and others, and it has become a starting 
point of many sophisticated approaches to fcrmionic 
many-body problems. 

For periodic systems, the effective free-electron the- 
ory and the band theory of Bloch [57] were the first 
crucial steps towards our present understanding of 
the real crystals. In 1930s, Wigner and Seitz [28l [29] 
performed the first quantitative calculations of the elec- 
tronic states in sodium metal. Building upon the homo- 
geneous electron gas model, the density-functional the- 
ory (DFT) was invented by Hohenberg and Kohn [30 
and further developed by Kohn and Sham 31j who 
formulated the local density approximation (LDA) for 
the exchange-correlation functional. These ideas were 
later elaborated by including spin polarization [32j , by 
constructing the generalized gradient approximation 
(GGA) [33 [31], and by designing a variety of orbital- 
dependent exchange-correlation functionals [35H37] . 
The DFT has proved to be very successful and has 
become the mainstream computational method for 
many applications, which cover not only solids but 
also molecules and even nuclear and other systems 
[38] 139] . The density- functional theory together with 
the Hartree-Fock and post-Hartree-Fock methods [40] 



are relevant for our discussion of quantum Monte Carlo 
methodology, since the latter uses the results of these 
approaches as a reference and also for construction 
of the many-body wave functions. Familiarity with 
the basic concepts of the Hartree-Fock and density- 
functional theories is likely to make the subsequent 
sections easier to follow, but we believe that it is not 
a necessary prerequisite for understanding our exposi- 
tion of the QMC methods and their foundations. 

2 Methods 

2.1 Variational Monte Carlo 

In the variational Monte Carlo method, the ground 
state of a hamiltonian H is approximated by some trial 
wave function I^t), whose form is chosen following a 
prior analysis of the physical system being investigated. 
Functional forms relevant to solid-state applications 
will be discussed later in section[4] Typically a number 
of parameters are introduced into I^t), and these pa- 
rameters are varied to minimize the expectation value 
£^2 = (* T |i?|4' T )/(* T |* T ) in order to bring the 
trial wave function as close as possible to the actual 
ground state \^o). 

Wave functions of interacting systems are non- 
separable, and the integration needed to evaluate £^2 
is therefore a difficult task. Although it is possible to 
write these wave functions as linear combinations of 
separable terms, this tactic is viable only for a limited 
number of particles, since the length of such expansions 
grows very quickly as the system size increases. The 
variational Monte Carlo method employs a stochastic 
integration that can treat the non-separable wave func- 
tions directly. The expectation value E^ is written 
as 



\V T (K)\ 2 [H*t](K) 



E 



1 



d 3N K 



VMC — 



A/" 

2— 1 



[h*t] {Hi) 



(3) 



where 1Z = (ri,V2, ■ ■ . , r/v) is a 3./V-dimensional vec- 
tor encompassing the coordinates of all N particles 
in the system and the sum runs over Af such vec- 
tors {7Zi} sampled from the multivariate probability 
density p(K) = |* t (^)| 2 /(*t|*t)- The summand 
E L {K) = [£f# T ](ft)/*T(ft) is usually referred to as 
the local energy. We assume spin-independent hamil- 
tonians, and therefore spin variables do not explic- 
itly enter the evaluation of the expectation value pi . 
This statement is further corroborated in section \A. II 
where the elementary properties of the trial wave func- 
tions |\Pt) are discussed. 
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Equation ^ transforms the multidimensional in- 
tegration into a problem of sampling a complicated 
probability distribution. The samples {TZi} can be 
obtained such that they constitute a Markov chain 
with transitions IZi+i <— TZi governed by a stochastic 
matrix M(lZ i+ i IZ4) whose stationary distribution 
coincides with the desired probability density p(TZ), 

p(K') = ( M{K' <- n)p{Tl) d 3N TZ for all TZ' . (4) 



After a period of equilibration, the members of the 
Markov sequence sample the stationary distribution 
regardless of the starting point of the chain, provided 
the matrix M(1Z' <— TZ) is ergodic. Inspired by the 
way the samples explore the configuration space, they 
are often referred to as walkers. 

The Markov chain can be conveniently constructed 
with the aid of the Metropolis method [31 BT] . 
The transition matrix is factorized into two parts, 
M(W 4r- Tli) = T(W 4r- Ki)A(K' <r- Ki), which cor- 
respond to two consecutive stochastic processes: A 
candidate TZ' for (i + l)-th sample is proposed accord- 
ing to the probability T(1Z' <— TZi) and this move is 
then either accepted with the probability A(1Z' <— TZi) 
or rejected with the probability 1 — A(7Z' <— TZi). If 
the move is accepted, the new member of the sequence 
is 7Z i+ i = 7Z' , otherwise it is !Z i+ i = TZi. The length 
of the chain is thus incremented in either case. The 
acceptance probability A(7Z' <— TZi), complementing 
some given T(7Z' TZi) and p(TZ) such that the sta- 
tionarity condition Q is fulfilled, is not unique. The 
choice corresponding to the Metropolis algorithm reads 



A(TZ' ^TZ)= min 



T(TZ <- TZ') p(TZ') 



' T(TZ' <- TZ) p(TZ) 



(5) 



and depends only on ratios of T and p. Consequently, 
normalization of the trial wave function I^t) is com- 
pletely irrelevant for the Monte Carlo evaluation of the 
quantum-mechanical expectation values. The freedom 
to choose the proposal probability T(TZ' <— TZi) can 
be exploited to improve ergodicity of the sampling, for 
instance, to make it easier to overcome a barrier of 
low probability density p separating two high-density 
regions. A generic choice for T(TZ' <— TZi) is a Gaus- 
sian distribution centered at TZi with its width tuned 
to optimize the efficiency of the sampling. 

The variational energy £Vmc is a stochastic vari- 
able, and an appropriate characterization of the ran- 
dom error £Vmc ~ is thus an integral part of the 
variational Monte Carlo method. When the sampled 
local energies E^fTZi) are sufficiently well behaved [32] , 
this error can be represented by the variance of i?vMC- 
In such cases, the error scales as Af^ 1 / 2 and is pro- 
portional to fluctuations of the local energy. Reliable 



estimation of the variance of £Vmc is a non-trivial 
affair since the random samples {TZi} generated by 
means of the Markov chain are correlated. These cor- 
relations are not known a priori and depend on the 
particular form of the transition matrix M that varies 
from case to case. Nevertheless, it is possible to esti- 
mate the variance without detailed knowledge of the 
correlation properties of the chain with the aid of the 
so-called blocking method [15] . 

The fluctuations of the local energy E^ are reduced 
as the trial wave function I'I't) approaches an eigen- 
state of the hamiltonian, and E^ becomes a constant 
when I'I't) is an eigenstate. In particular, it is crucial 
to remove as many singularities from E^ as possible 
by a proper choice of the trial function. Section |4T] il- 
lustrates how it is achieved in the case of the Coulomb 
potential that is singular at particle coincidences. 

The total energy is not the only quantity of inter- 
est and evaluation of other ground-state expectation 
values is often desired. The formalism sketched so far 
remains unchanged, only the local energy is replaced 
by a local quantity A L (TZ) = [A* T ] (TZ)/^ T (TZ) cor- 
responding to a general operator A. An important 
difference between A^ and the local energy is that 
fluctuations of A^ do not vanish when I^t) is an 
eigenstate of H. These fluctuations can severely im- 
pact the efficiency of the Monte Carlo integration in 
(5't|^1| x E , t)/(^'t|*t), and the random error can decay 
even slower than A/" -1 / 2 [12]. The trial wave function 
cannot be altered to suppress the fluctuations in this 
case, but a modified operator A' can often be con- 
structed such that (* T |A'|* T ) = (*t|^|*t) while 
the fluctuations of A^ are substantially reduced [44T - 



2.2 Diffusion Monte Carlo 

The accuracy of the variational Monte Carlo method is 
limited by the quality of the trial wave function |\&t)- 
This limitation can be overcome with the aid of the 
projector methods. In particular, the diffusion Monte 
Carlo method [HI employs an imaginary time 

evolution 



|tf D (<)> =exp(-[#-Sr(i)]t) |* n 



(6) 



where the energy offset Et is introduced to main- 
tain the wave-function norm at a fixed value. Formal 
properties of (|6| can be elucidated by expanding the 
trial function |^t) m terms of the hamiltonian eigen- 
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states ([2]), which readily yields 
|* D (t)) = exp(-[E - E T (t)]t 



|*o)(*o|*r> 



J2e-^- E «^ n )(^ n \^ T ) 



(J) 



The ground state |^o) is indeed reached in the limit of 
large t as long as the trial function was not orthogonal 
to |^o) from the beginning. The requirement of a 
finite norm of (*)) translates into a formula 



E = lim E T (t) 

t—too 



(8) 



that can be used to obtain the ground-state energy. An 
alternative approach is to evaluate the matrix element 
£* D * T = (*dW|£|*t)/(*d(*)I*t) that asymptot- 
ically coincides with the ground-state energy, since 

(* |ff|*r)/(*o|*T> = (*o|£|*o)/(*o|*o)- The in- 
tegration in Eq, D \z, T can be performed stochastically 
in analogy with the VMC method, 



E 



(*d(*)|*t> 
~ Ebmc 



*t(^) 



1 N 

1—1 



E h (Ki), (9) 



where the individual samples IZi are now drawn 
from a probability distribution defined as p(1Z, t) = 

*D^iWR)/( $ D(t)|*T). 

2.2.1 Fixed-node/ fixed-phase approximation 

The Monte Carlo integration indicated in Q is pos- 
sible only if p(!Z, t) is real- valued and positive. Since 
the hamiltonians we usually deal with are symmetric 
with respect to time reversal, the eigenfunctions can 
be chosen real. Unfortunately, many-electron wave 
functions must necessarily have alternating sign to 
comply with the fermionic antisymmetry. In general, 
the initial guess |\I/t) will have different plus and mi- 
nus sign domains (also referred to as nodal pockets 
or nodal cells) than the sought for ground-state wave 
function |\&o)j which results in changing sign of p(1Z, t). 
In certain special cases, the correct sign structure of 
the ground state can be deduced from symmetry con- 
siderations [52H54] , but in a general interacting system 
the exact position of the boundary between the pos- 
itive and negative domains (the so-called fermionic 
node) is unknown and is determined by the quantum 
many-body physics |55j . A number of exact properties 
of the fermionic nodes have been discovered [56l - f59j . 
but a lot remains to be done in order to transform this 



knowledge into constructive algorithms for the trial 
wave functions. 

The problem with the variable sign of p(TZ, t) can 
be circumvented by complementing the projection ([6]) 
with the so-called fixed-node constraint [15] , 

t)V T {K) > for all K and all t . (10) 

Doing so, lim^oo |\£d(£)) only approximates |^o), 
since the projection cannot entirely reach the ground 
state if the initial wave function \^t) does not possess 
the exact nodes. The total energy calculated with 
this fixed-node method represents an upper-bound es- 
timate of the true ground-state energy because the 
projection ^ is restricted to a subspace of the whole 
Hilbert space when the constraint ( 10 ) is implemented 



The fixed-node approximation has proved 
very fruitful in quantum chemistry |631 164) as well as 
for investigation of the electronic structure of solids 
as testified by the applications reviewed in section [5j 
In calculations of extended systems and espe- 
cially metals, it is beneficial to allow for boundary 
conditions that break the time-reversal symmetry, 
since it facilitates reduction of finite-size effects (sec- 
The eigenfunctions are then complex- valued 



3.1) 



tion 

and a generalization of the fixed-node approxima- 
tion is required 



The constraint ( 10 ) is replaced with 
vp D (i) = \^ u (t)\e iipT , where ip T is the phase of the 
trial wave function ^>t = \^t \ e 1VT [65] ■ The phase <^t 
is held constant during the DMC simulation to guar- 
antee that p(lZ, t) stays non-negative for all TZ and t. 
Additionally, a complex trial wave function |^t) causes 
the local energy E^ to be complex as well. The ap- 
propriate modification of the estimate for the total en- 
ergy (|9| coinciding with the asymptotic value of E^{t) 
then reads 

1 M 

£ D MC = ^^Re[£ L (7^)]. (11) 



Analogous to the fixed-node approximation, the fixed- 
phase method provides a variational upper-bound es- 
timate of the true ground-state energy. Moreover, 
the fixed-phase approximation reduces to the fixed- 
node approximation when applied to real-valued wave 
functions. 

2.2.2 Sampling the probability distribution 

The unnormalized probability distribution that we 
wish to sample in the fixed-phase DMC method, 



f(K,t) = $i(K,i)$ T (R) = |*D(ft,t)||*T(ft)| 



(12) 



G 



referred to as the mixed distribution, fulfills an equa- 
tion of motion 



dtf(K,t) = 

-^V 2 f(K,t)+V-[v B (K)f(K,t)] 

+ f(TZ,t)\Re[E L (TZ)] - (1 + td t )E T {t)] (13) 



that is derived by differentiating Q and (12) with 
respect to time, combining the resulting expressions 
and rearranging the terms. The drift velocity in- 
troduced in (13) is defined as i>d = V1h|$t| an d 



V denotes the 3-ZV-dimensional gradient with respect 
to 1Z. The equation of motion is valid in this form 
only as long as the kinetic energy is the sole non-local 
operator in the hamiltonian. Strategies for inclusion 
of non-local pseudopotentials will be discussed later 
in section |2.3| The case of the fixed-node approxi- 



mation is virtually identical to (13), except that the 



local energy is real by itself. The following discussion 
therefore applies to both methods. 

The time evolution of the mixed distribution 
f(JZ,t) can be written in the form of a convolution 



f(TZ,t) = J 'g(K<-K' ,t)f(W \0)d 3N TZ / . 



(14) 



where f(lZ,0) = | \I/ t (^-) 1 2 and the Green's function 
G{K <- 72/, i) = (Tl\G{t)\K') is a solution of @ 
with the initial condition G(K <- W , 0) = 5(K - K'). 
Making use of the Trotter-Suzuki formula [5^1 Wf\ , 
the Green's function is approximated by a product of 
short-time expressions, 

G{t) = [G g / d (r) G diff (r) G drift (r)] M + O(r) , (15) 

" v ' 

G st (r) 



where r denotes t/M and the exact solution of (13) 
is approached as this time step goes to zero. Con- 
sequently, the DMC simulations should be repeated 
for several sizes of the time step and an extrapola- 
tion of the results to r — > should be performed in 
the end. For simplicity, we show in (15) only the 



simplest Trotter-Suzuki decomposition which has a 
time step error proportional to r. Commonly used are 
higher order approximations whose errors scale as r 2 
or t 3 . The three new Green's functions constituting 
the short-time approximation G s t can be explicitly 



written as 

G dim (Tl^n',r) (16) 
= [1 - rV • v u (K')]8[K - W - v d (K')t] 
+ 0(r 2 ), 



G diff {n<- n', t) 
i 

- (27tt) 3A 7 2 
G g/d (TZ^n',T) 



(17) 



exp 



(K - 1l') 2 
2t 



(18) 



-T[Re[E L (n)] - E T (ty) s[n-n'] 



— exp 

and correspond to the three non-commuting oper- 
ators from the right-hand side of ( [13] ) in the or- 
der: drift (V • [v u (U)»]), diffusion (^V 2 /2.) and 
growth/decay (•[Re[E L (7Z)] - (1 + td t )E T (t)]). The 
drift and diffusion Green's functions preserve the nor- 
malization of f(lZ,t) whereas the growth/decay pro- 
cess does not. 

The factorization of the exact Green's function into 
the product of the short-time terms forms the basis 
of the stochastic process that represents the diffusion 
Monte Carlo algorithm. First, A4 samples {IZi} are 
drawn from the distribution f(K,0) = |*t(^)| 2 just 
like in the VMC method. Subsequently, this set of 
walkers evolves such that it samples the mixed distri- 
bution f(JZ,t) at any later time t. The probability 
distribution is updated from time t to t + t by multi- 
plication with the short-time Green's function, 

f(K,t+r)= f G st (R^K>,T)f(K',t)d 3N K' , (19) 



which translates into the following procedure per- 
formed on each walker in the population: 

(i) a drift move AT^drift = vd(1Z')t is proposed 

(ii) a diffusion move ATZ^g — x is proposed, where 
X is a vector of Gaussian random numbers with 
variance r and zero mean 



(iii) 



the increment A72-drift + AT^diff is accepted 
if it complies with the fixed-node condition 
*t(^')*t(^' + Adrift + Aft di ff) > 0, other- 
wise the walker stays at its original position; 
moves attempting to cross the node occur only 
due to inaccuracy of the approximate Green's 



function ( 15 ), and they vanish in the limit t — > 0: 
the moves A7\Ld r ift + A7\Ldiff are accepted without 
any constraint in the fixed-phase method 

(iv) the growth/decay Green's function G g / d is ap- 
plied; several algorithms devised for this purpose 
arc outlined in the following paragraph 
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(v) at this moment, the time step is finished and the 
simulation continues with another cycle starting 
back at 

After the projection period is completed, the algorithm 
samples the desired ground-state mixed distribution 
and the quantities needed for evaluation of various 
expectation values can be collected in step Q. 

At this point we return to a more detailed 
discussion of several algorithmic representations of 
the growth/decay Green's function G g / d needed in 
step ( pv| . 

• The most straightforward way is to assign 
a weight w to each walker. These weights 
are set to 1 during the VMC initialization 
of the walker population and the applica- 
tion of Gg/d then amounts to a multiplication 
w(t + t) = w(t)W(1Z), where the weighting fac- 
tor is 

W{K) = exp -T(Re[E h (K)] -E T (t)) . (20) 
Consequently, the formula for calculation of the 



total energy (11) is modified to 



E 



DMC 



Wj 



-l M 

£ 

i=l 



Re[E L (Ki)] (21) 



and the walkers remain distributed according 
to \^ T (TZ)\ 2 as in the VMC method. This algo- 
rithm is referred to as the pure diffusion Monte 
Carlo method [68 , 69 . It is known to be intrinsi- 
cally unstable at large projection times where the 
signal-to-noise ratio deteriorates [70], but it is 
still useful for small quantum-chemical systems 

[SHS|. 

• The standard DMC algorithm fixes the weights 
to w = 1 and instead allows for stochastically 
fluctuating size of the walker population by 
branching walkers in regions with large weight- 
ing factor W(1Z) and by removing them from 
areas with small W(7Z). The copies from high- 
probability regions are treated as independent 
samples in the subsequent time steps. The time 
dependence of the energy offset Er(t) provides a 
population control mechanism that prevents the 
population from exploding or collapsing entirely 
[501 [73]. The branching/elimination algorithm is 
much more efficient in large many-body systems 
than the pure DMC method, although it also 
eventually reaches the limits of its applicability 
for a very large number of particles [75) . 



• An alternative to the fluctuating population are 
various flavours of the stochastic reconfigura- 
tion [T5ir7Dl l75H7gl. These algorithms comple- 
ment each branched walker with high weighting 
factor W(1Z) with one eliminated walker with 
small W(TZ), and therefore the total number of 
walkers remains constant. This pairing intro- 
duces slight correlations into the walker popu- 
lation that are comparable to those caused by 
the population control feedback of the standard 
branching/elimination algorithm [75] , Keeping 
the population size fixed is advantageous for 
load balancing in parallel computational envi- 
ronments, since the number of walkers can be 
a multiple of the number of computer nodes 
(CPUs) at all times during the simulation. 



The branching/elimination process interacts in a subtle 
way with the fixed- node constraint. Since the walk- 
ers are not allowed to cross the node, the branched 
and parent walkers always stay in the same nodal cell. 
If some of these cells are more favoured (that is, if 
they have a lower local energy on average), the walker 
population accumulates there and eventually vanishes 
from the less favoured cells. Such uneven distribution 
of samples would introduce a bias to the simulation. 
Fortunately, it does not happen, since all nodal cells 
of the ground-state wave functions are connected by 
particle permutations and are therefore equivalent, see 
the tiling theorem in [SO] . For general excited states 
this theorem does not hold and the unwanted depopu- 
lation of some nodal cells can indeed be observed. The 
problem is absent from the fixed-phase method, since 
it contains no restriction on the walker propagation. 

The branching/elimination algorithm is just one 
of the options of dealing with the weights along the 
stochastic paths. Another possibility was introduced 
by Baroni and Moroni as the so-called reptation al- 
gorithm [79 , which recasts the sampling of both the 
path in the configuration space and the weight into 
a straightforward Monte Carlo process, avoiding thus 
some of the disadvantages of the DMC algorithm. The 
reptation method has its own sources of inefficiencies 
which can be, however, significantly alleviated [BO] . 

This concludes our presentation of the stochastic 
techniques that are used to simulate the projection 
operator introduced in Q. We would like to bring 
to the reader's attention that the algorithm outlined 
in this section is rather rudimentary and illustrates 
only the general ideas. A number of important perfor- 
mance improvements are usually employed in practical 
simulations, sec for instance [74] for further details. 
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2.2.3 General expectation values 

So far, only the total energy was discussed in con- 
nection with the DMC method. An expression anal- 
ogous to ([9]) can be written with any operator A 
in place of the hamiltonian H. The acquired quan- 
tity A^ u q, T = (* d |A|* t )/(*d| , I't), called the mixed 
estimate, differs from the pure expectation value 
(\]/d|-4| v I / d)/( v I'd| v I / d> unless A commutes with the 
hamiltonian. In general, the error is proportional to 
the difference between \^d) an d \*&t)- The bias can 
be reduced to the next order using the following ex- 
trapolation [71150] 



(*D\A\*n) 

(#£>|#d) 



(*_d|^4|*t) {9t\A\^t) 



O 



(*d|*t) 



W T W T ) 



(22) 



Alternative methods that allow for a direct evaluation 
of the pure expectation values have been developed, 
such as the forward (or future) walking J50j [HJ HU , 
the reptation quantum Monte Carlo [7j2 [83j [84] , or the 
Hellman-Feynman operator sampling [85, 86 . Due 
to their certain limitations, these techniques do not 
fully replace the extrapolation ( 22 ) — they are usable 



only for local operators and the former two become 
computationally inefficient in large systems. 

The discussion of the random errors from the end 
of section [2~T] applies also to the diffusion Monte Carlo 
method, except that the serial correlations among the 
data produced in the successive steps of the DMC sim- 
ulations are typically larger than the correlations in 
the VMC data. Therefore, longer DMC runs are neces- 
sary to achieve equivalent suppression of the stochastic 
uncertainties of the calculated expectation values. 

2.2.4 Spin degrees of freedom 

The DMC method as outlined above samples only the 
spatial part of the wave function, and the spin degrees 
of freedom remain fixed during the whole simulation. 
This simplification follows from the assumption of a 
spin-independent hamiltonian that implies freezing of 
spins during the DMC projection ([6|. This is indeed 
the current state of the DMC methodology as applied 
to electronic-structure problems: In order to arrive at 
the correct spin state, a number of spin-restricted cal- 
culations are performed and the variational principle 
is employed to select the best ground state candidate 
among them. 

Fixing spin variables is not possible for spin- 
dependent hamiltonians, such as for those containing 



spin-orbital interactions, since they lead to a non- 
trivial coupling of different spin configurations. In 
fact, spin-dependent quantum Monte Carlo methods 
were developed for studies of nuclear matter. A variant 
of the Green's function Monte Carlo method [87, 88 
treats the spin degrees of freedom directly in their full 
state space. Since the number of spin configurations 
grows exponentially with the number of particles, this 
approach is limited to relatively small systems. More 
favourable scaling with the systems size offers the aux- 
iliary field diffusion Monte Carlo method that samples 
the spin variables stochastically by means of auxiliary 
fields introduced via the Hubbard-Stratonovich trans- 
formation (SHI HO] • Recently, a version of the auxiliary 
field DMC method was used to investigate properties 
of the two-dimensional electron gas in presence of the 
Rashba spin-orbital coupling [9"T] . 

2.3 Pseudopotentials 

The computational cost of all-electron QMC calcula- 
tions grows very rapidly with the atomic number Z 
of the elements constituting the simulated system. 
Theoretical analysis [SSI [S3] as well as practical cal- 
culations [53] indicate that the cost scales as z 5 - 5 ~ 6 5 . 
Most of the computer time spent in these simulations 
is used for sampling of large energy fluctuations in the 
core region, which have very little effect on typical 
properties of interest, such as interatomic bonding 
and low-energy excitations. For investigations of these 
quantities it is convenient to replace the core electrons 
with accurate pseudopotentials. A sizeable library of 
norm-conserving pseudopotentials targeted specifically 
to applications of the QMC methods has been built 
over the years [9"5TU00| . 

Pseudopotentials substitute the ionic Coulomb po- 
tential with a modified expression, 



Z 
r 



V(r) + W 



(23) 



where V(r) is a local term behaving asymptotically as 
— (Z — Z covc )/r with Z corc being the number of elim- 
inated core electrons. The operator W is non-local 
in the sense that electrons with different angular mo- 
menta experience different radial potentials. Explicitly, 
the matrix elements of the potential W associated with 
7-th atom in the system are 



(n\w I \n') = J2z2 



1 



i 1—0 m=—l 



[fu\lm) 



x Wmiri^im-r'^ilrrAr'ii), (24) 

where \lm) are angular momentum eigenstates, rn is 
the distance of an electron from the 7-th nucleus and 
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rn is the associated direction Tu/ru. Functions Wjj 
vanish for distances ru larger than some cut-off ra- 
dius r c , and the sum therefore runs only over 
electrons that are sufficiently close to the particular 
nucleus. 

Evaluation of pseudopotentials in the VMC method 
is straightforward, despite the fact that the local en- 
ergy Ei, itself involves integrals. As can be inferred 



from the form of the matrix elements (24 1, these are 



two-dimensional integrals over surfaces of spheres cen- 
tered at the nuclei. The integration can be imple- 
mented with the aid of the Gaussian quadrature rules 
that employ favourably sparse meshes [1011 1102J. 

The use of non-local pseudopotentials in the fixed- 
node DMC method is more involved, since the sam- 
pling algorithm outlined in section [2 . 2 . 2| explicitly as- 
sumed that all potentials were local. Non-local hamil- 
tonian terms can be formally incorporated by introduc- 



ing an extra member into the Trotter break-up (15), 
namely 

G nloc (K^K',T) 

= W - - (n\rW\W) + 0(r 2 ) , (25) 

where W now combines the non-local contributions 
from all atoms in the system. This alone is not the 
desired solution, since the term involving the matrix 
element of W does not have a fixed sign and thus 
cannot be interpreted as a transition probability. 

To circumvent this difficulty the so-called localiza- 
tion approximation has been proposed. It amounts to 
a replacement of the non-local operator in the hamil- 
tonian with a local expression [931 11021 1103] 

* - ^-W- <26 » 

Consequently, the contributions from W are directly 
incorporated into the growth/decay Green's func- 
tion ( 18 1 and no complications with alternating sign 



arise. Unfortunately, the DMC method with this ap- 
proximation does not necessarily provide an upper- 
bound estimate for the ground-state energy. The cal- 
culated total energy Eumc is above the lowest eigen- 
value of the localized hamiltonian, which does not 
guarantee that it is also higher than the ground-state 
energy of the original hamiltonian H. The errors in 
the total energy incurred by the localization approx- 
imation are quadratic in the difference between the 
trial function I^t) an d the exact ground-state wave 
function [102] . The trial wave functions are usually 



accurate enough for the localization error to be practi- 
cally insignificant and nearly all applications listed in 
section [5] utilize this approximation. 

A method that preserves the upper-bound property 
of -Edmc was proposed in the context of the DMC algo- 
rithm developed for lattice models [104] . The non-local 
operator W is split into two parts, W = W + + W~ , 
such that W + contains those matrix elements, for 
which (ft|W"|72, / )*T(ft)*T(ft / ) is positive, and W~ 
contains the elements, for which the expression is neg- 
ative. Only the W + part is localized in order to obtain 
the approximate hamiltonian, 



W 



w- 



$ T (K) 



(27) 



One can explicitly show that the lowest eigenvalue of 
this partially localized hamiltonian is an upper bound 
to the lowest eigenvalue of the original fully non-local 
hamiltonian [104 . Recently, a stochastic represen- 



tation of the non-local Green's function (25) corre- 
sponding to W~ was implemented also into the DMC 
method for continuous space [105J. Apart from the 
recovered upper-bound property, the new algorithm 
reduces fluctuations of the local energy for certain 
types of pseudopotentials. On the other hand, the 
time step error is in general larger [1051 1106j , since the 
distinct treatment of the W + and W~ parts of the 
pseudopotential essentially corresponds to a Trotter 
splitting of the growth/decay Green's function (18) 
into two pieces. Very recently, a more accurate Trotter 
break-up and other modifications improving efficiency 
of this method have been proposed for both continuous 
and lattice DMC formulations [107] , 

The localization approximation is directly applica- 
ble also to the fixed-phase DMC method. Adaptation 
of the non-local moves representing W~ to cases in- 
volving complex wave functions has not been reported 
yet, nevertheless, the modifications required should be 
only minor. 



3 From a finite supercell to the 
thermodynamic limit 

Quantum Monte Carlo methods introduced in the pre- 
ceding chapter can be straightforwardly applied to 
physical systems of a finite size, such as atoms and 
clusters of atoms. To allow investigation of bulk prop- 
erties of solids, the algorithms described so far have to 
be complemented with additional techniques that re- 
duce the essentially infinitely many degrees of freedom 
into a problem of manageable proportions. 
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3.1 Twist-averaged boundary conditions 



and of an electron-electron contribution 



In approximations that model electrons in solids as an 
ensemble of independent (quasi-)particles, it is possible 
to map the whole infinite crystal onto a finite volume 
so that the thermodynamic limit becomes directly ac- 
cessible. Hamiltonians of such models are invariant 
with respect to separate translations of electrons by 
any lattice vector _R, that is, for each i we can write 



H{r 1 ,r 2 , 



R, 



H(r 1 ,r 2 , 



.). (28) 



This invariance allows to diagonalize the hamiltonian 
only in the primitive cell of the lattice and then use 
the translations to expand the eigenstates from there 
into the whole crystal. Unfortunately, the explicit 
two-body interactions that we are set out to keep in 



the hamiltonian break the symmetry ( 28 1 . The only 



translation left is a simultaneous displacement of all 
electrons by a lattice vector, which is not enough to 
reach the thermodynamic limit with a finite set of 
degrees of freedom. 

To proceed further we introduce artificial trans- 
lational symmetries with the aid of the so-called 
supercell approximation that is widely used within 
the independent-particle methods to investigate non- 
periodic structures such as lattice defects. We select 
a supercell having a volume Q$ that contains several 
(preferably many) primitive cells. The whole crystal 
is then reconstructed via translations of this large 
cell by supercell lattice vectors {-Rs}, which are a 
subset of all lattice vectors {R}. Simultaneously, we 
divide the electrons in the solid into groups containing 
N = p av Qs particles, where p av is the average electron 
density in the crystal. This partitioning is used to 
construct a model hamiltonian, where electrons within 
each group interact, whereas there are no interactions 
between the groups, 



i=i 



E 

i=i 



N 



EW+^.(w (1) ) 



(29) 



The vector = (rj , . . . , rjy ) denotes coordinates 
of the electrons belonging to the 7-th group. Note 
that these electrons are not confined to any particular 
region in the crystal. The supercell hamiltonian H$ 
consists of single-particle terms h, which encompass ki- 
netic energy as well as ionic and all external potentials, 



E 



1 E 

2 ^ 



Rs 



(30) 



The first term in ( 30 ) represents the explicit Coulomb 



interaction among electrons in the iV-member group, 
and the second term mimics the interactions with the 
electrons outside the group. The physical meaning 
of the latter term is as follows: A set of images is 
associated with each electron j, and these virtual par- 
ticles are placed at positions rj + Rs so that they 
create a regular lattice. The combination of all images 
has the same average charge density as the original 
crystal and thus represents a reasonable environment 
for the selected N electrons. Each electron i then 
interacts with the arrays of charges associated with 
the other electrons in the group as well as with its 
own images. Only one half of the interaction energy 



with images is included in ( 30 1 , the other half belongs 



to the rest of the system and is distributed among the 



other members of the sum in ( 29 ) . The model hamilto- 
nian H m odci approaches the original fully interacting 
hamiltonian as N increases and a larger fraction of 
the interactions has the exact form. 

Hamiltonians Hg and i? m odci possess the symmetry 



described by ( 28 ) with lattice translations R replaced 



with Rs- Consequently, the eigenfunctions of Hs are 
many-particle Bloch waves 



*K a (R) = U Ka (TZ)exp [iK 



N 



(31) 



where a is a many-body analogue of the band in- 
dex and K is the crystal momentum |1081 1109] . The 
wave functions of the form (31 1 can be found in the 



same way as the single-particle Bloch waves within 
the independent-particle methods — as solutions to a 
problem of N particles confined to a simulation cell 
defined by vectors Li, L 2 and L 3 belonging to the 
set {Rs} and giving fig = \Li ■ (L 2 X L 3 )\. The dy- 
namics of this finite ./V-particle system is governed by 
the hamiltonian Hs complemented with the so-called 
twisted boundary conditions |110j 

®K a {ri, ■ • ■ ,n + L m , . . .,r N ) 

= *K a (ri,...,ri,...,r N )e iK - L ™ , (32) 

where i = l,...,N and m = 1,2,3. The indistin- 
guishability of electrons implies that the phase factor 
is the same irrespective of which electron is moved, 
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Figure 1: Deviation of the twist-averaged to- 
tal energy ( 35 I from the exact thermodynamic 
limit E x , A = [Es(N) - E^/N, for a three- 
dimensional gas of non-interacting electrons with 
dispersion relation = fc 2 /2 at density p corre- 
sponding to r s = [3/(47rp)] 1/3 = 1. The dashed 
line represents the average asymptotic decay of 
A that behaves as iV -1 ' 32 . 



> 




N 



and therefore only a single K vector common to all 



electron coordinates is allowed in (31 1 and (32). Once 



all twisted boundary conditions [11011111] . In practice, 
the integral in ( 35 ) is approximated by a discrete sum. 



the quantum-mechanical problem in the finite simula- 
tion cell is solved, wave functions for the whole crystal 
can be constructed. Since there are no interactions 
between individual iV-particle groups, these wave func- 
tions have the form of an antisymmetrized product of 



The larger the simulation cell the smaller number of 
K points is needed to represent the integral, since the 
first Brillouin zone of the simulation cell shrinks and 
the ^-dependence of the integrand gets weaker with 
increasing fig. 



the Bloch functions (31), namely 



Formula ( 35 ) is almost identical to the expression 



{K l}{ai }{{r}) =A 



n 

i=i 



Krai 



(33) 



The indicated antisymmetrization is straightforward 
as long as all Kj in the product are different, because 
each factor $jf faj then comes from a disjoint part 
of the Fock space. The total energy corresponding 



used in supercell calculations within the independent- 
particle theories, the only difference is that the number 
of electrons at each K point is fixed to N. This con- 
straint is benign in the case of insulators where the 
number of occupied bands is indeed constant across 
the Brillouin zone. In metals, however, the number 
of occupied bands fluctuates from one K point to an- 



to the wave function ( 33 1 with the extra restriction 
Kj Kp reads 

oo 

E {Kl7 t Kl ,}{ ai } = ^2{9/K iai \H s \^K iai ) (34) 
i=i 

and the lowest energy is obtained by setting aj = 0, 
that is, by selecting the ground state for each of the 
different boundary conditions. Although unlikely, it 
is possible that the true ground state of H moAc \ falls 
outside the constraint Ki ^ Kp. In those cases, the 



other, and therefore the average (35 1 does not coincide 



E S = (H S ) 



(2tt) 2 



l.B.Z. 



with the exact thermodynamic limit. Figure [T] pro- 
vides an illustration of the residual error. In principle, 
it is possible to remove this error with the aid of the 
grand-canonical description of the simulation cell |110j , 
but this concept is not straightforward to apply since 
the supercell is no longer charge neutral. 

3.2 Ewald formula 

Our definition of the simulation-cell hamiltonian Hg 
in the preceding section was only formal and deserves 
a further commentary. It turns out that V ee is not 
absolutely convergent, and therefore it does not un- 
ambiguously specify the interaction energy. In partic- 
ular, the seemingly periodic form of the sums in (30) 
does not by itself imply the desired periodicity of the 
hamiltonian. However, enforcing this periodicity as 
an additional constraint makes the definition unique 
and the resulting quantity is known as the Ewald en- 

d 3 K (^ ko\Hs\^ ko) (35) er Sy ^ e ■ ^ can ^ e s h° wn that the requirement of 

periodicity is equivalent to a particular boundary con- 
dition imposed on the electrostatic potential at infinity 
|112( 1113] . The peculiar convergence properties of the 
sums in ( 30 1 are a manifestation of the long-range 



expression (34) with on — is an upper-bound esti- 
mate of the actual ground-state energy of -ff mo dci with 
a bias diminishing as iV increases. Taking into account 
the continuous character of the momentum K in the 
infinite crystal and the fact that all possible boundary 



conditions ( 32 ) are exhausted by all K vectors within 



the first Brillouin zone, the ground-state energy per 
simulation cell can be written as 



The total energy as well as expectation values of other 
periodic operators are calculated as an average over 
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character of the Coulomb potential — the boundary of 
the sample is never irrelevant, no matter how large 
its volume is. Consequently, careful considerations are 
required in order to perform the thermodynamic limit 
correctly. 

For the purposes of practical evaluation in QMC 
simulations, the Ewald energy is written as 



an appropriate function e fit (iV) = + f(N) is subse- 
quently fitted through the acquired data. In the end, 
Eoo is the desired energy per electron in the thermody- 
namic limit, where the error term f(N) by definition 
vanishes, limjv->.oo f(N) = 0. Experience with a wide 
range of systems [TU1 I117H119] indicates that as long 
as the integral over the Brillouin zone in ( 35 ) is well 



K ( e E) W = E^(n-r,) 

j<i 

+ -lim(y E (r) 



(36) 



where Ve(^« — Tj) stands for an electrostatic potential 
at Tj induced by the charge at rj together with its 
images located at rj + Rs- An explicit formula for 
the Ewald pair potential reads [1121 1114] 



V E (ri-rj) 



Rs 1 
4tt 



— — erfc(K|r l - r 3 - R s \) 



+ 



Gs^O 



converged, the finite-size data are well approximated 
by a smooth function f(N) < dominated by a 1/N 
contribution)^] The size extrapolation is therefore quite 
straightforward, although often computationally ex- 
pensive due to the relatively slow decay of the error 
term. 

The origin and behaviour of the finite-size effects 
is a subject of ongoing investigations with the aim 
to find means of accelerating the convergence of the 
total energy and other expectation values to the ther- 
modynamic limit. Furthermore, understanding the 
dependence of the finite-size errors on various parame- 
ters, such as particle density, can reduce the number of 
explicit size extrapolations needed to obtain quantities 
of interest. In calculations of equations of state (sec- 
tion 5.3), for instance, it is then sufficient to perform 



(37) 



where Gs are vectors reciprocal to Rs, exp(iGs ■ 
Rs) = 1, and k is an arbitrary positive constant that 
does not alter the value of Ve- The omission of the 
Gs = term in the reciprocal sum corresponds to 
the removal of the homogeneous component from the 
potential Ve- When evaluating the total energy of 
a charge neutral crystal, these "background" contri- 
butions exactly cancel among the Ewald energies of 
electron-electron, electron-ion and ion-ion interactions. 



the extrapolation only at selected few, instead of all, 
electron densities |119j . 

It turns out that, in the twist-averaged expecta- 
tion values (Hs) calculated in finite simulation cells, 
both the hamiltonian and the wave function contain 
biases that contribute to the 1/N asymptotics of the 
error term f(N). It was argued [1131 1120] that the 
slow converging parts of the hamiltonian reside in the 
exchange-correlation energy 



V 



xc 



p(r)VE(r-r')p{r')d 3 rd 3 r' (38) 



The important feature of the Ewald formula ((37} is defined as a difference between the expectation value 

of the Ewald energy (Vee } and the Hartree term 
that describes the interaction of the charge densities 
p{r) = (p(r)}. The Hartree energy is found to con- 
verge rather rapidly with the size of the simulation 
cell. In systems with cubic and higher symmetry, the 
leading contribution to f(N) can be written as |121| 



the decomposition of the slowly converging Coulomb 
sum into two rapidly converging parts, one in real 
space and the other in reciprocal space. The break- 
up is not unique (not only due to the arbitrariness 
of k) and can be further optimized for computational 
efficiency [TI51 Eg. 

3.3 Extrapolation to the thermodynamic 
limit 

The total energy per electron ejv = Es/N evaluated 
according to the outlined recipe still depends on the 
size of the simulation cell. These residual finite-size 
effects can be removed by an extrapolation: Energy ejv 
is calculated in simulation cells of several sizes and 



/xc(iV) 



Vxc ,. Vxc 

— — hm — — 

N N->oo N 

= _2tt ^ S(G S ) 
Q s Gs^o G 2 S 

This formula involves the static structure factor 



(39) 



S(G S ) = i [(p(Gs)p(-Gs)) - \(P(G S ))\ 



(40) 



1 Thc finite-size scaling depends on the dimensionality of the problem and the 1/N dependence corresponds to three-dimensional 
samples considered here. 
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where p(Gg) is a Fourier component of the density 
operator. The exact small-momentum asymptotics of 
the structure factor in Coulomb systems is given by the 
random phase approximation (RPA) |122j and reads 



S(Gs) ~ G|, which ensures that the limit in (39) is 



well defined. In systems with lowered symmetry and 
for less accurate approximations, as is the Hartree- 
Fock theory where S(Gs) ~ Gs, the expression for 
fxc(N) must be appropriately modified jl23) . 

The random phase approximation provides insight 
also into the finite-size effects induced by restricting 
the wave function into a finite simulation cell. Ac- 
cording to the RPA, the many-body wave function in 
Coulomb systems factorizes as 



V(n) = * B . r .(ft)exp 



(41) 



where 'J's.r. contains only short-range correlations and 
the function u(r) decays as 1/r at large distances. 
Such long-range behaviour is not consistent with the 



boundary conditions (32), and a truncation of this tail 



is therefore unavoidable. The corresponding finite-size 
bias is most pronounced in the expectation value of 
the kinetic energy T = (T) and contributes an error 
term [T2"T] 



T 

N \ 
1 



lim 



T 

N 



„ lim G 2 s u(G s ) 



(42) 



where u(Gs) ~ 1/G| is a Fourier component of u(r). 

Assuming that we are able to evaluate the expres- 
sions p9| and (42), we can decompose f(N) into 
parts f(N) = f xc (N) + f T (N) + f'(N), where f(N) 
is substantially smaller than f(N) and the size extrap- 
olation is therefore better controlled. In the case of 
the homogeneous electron gas, the RPA provides ana- 
lytic expressions for the small momentum behaviour 
of the required quantities S(Gs) ~ G|/(2w p ) and 
u(G s ) « in/(GgLJ p ), where uj p = ^AirN/^l s is the 
plasma frequency. Subsequently, the individual error 
terms simplify to 



(43) 



It can be demonstrated that these two contributions 
completely recover the 1/2V part of f(N), and that 
the residual term f(N) scales as ~ 2V" 4 / 3 [T2"3"] . 

Application of the derived finite-size corrections to 
simulations of realistic solids is less straightforward 
since reliable analytic results are not available. The 
small momentum asymptotics of S(Gs) and u(Gs) 



have to be examined numerically. Sufficiently large 
simulation cells are needed for this purpose, since the 
smallest nonzero reciprocal vector available in a su 
percell with volume fig is Gs 



the kinetic energy correction (42) within DMC sim 



flg 1 ^ 3 . Utilization of 



ulations is further complicated by the fact that the 
function u(r) is not given as an expectation value of 
an operator, and thus it is not clear how it could be 
extracted from the sampled mixed distribution \f r J ) \I , T- 
One has to rely on the trial wave function alone to 



correctly reproduce the long-range tail (41 ), which can 



be a challenging task in simulation cells containing a 
large number of electrons. 

The above analysis employs exact analytic formu- 
lae or quantum Monte Carlo simulations themselves 
to find corrections to the finite-size biases. Alterna- 
tively, it is possible to adopt a more heuristic approach 
and estimate the finite-size effects within an approx- 
imative method. For instance, the error term f(N) 
can be rewritten as f(N) = e s LDA) - e£ DA) + /"(2V), 
where e s LDA ' ) and 6^°^ are total energies per particle 
provided by the local density approximation with two 
distinct exchange-correlation functionals, and f"(N) is 
anticipated to be considerably smaller than f(N) |124j . 
The exchange-correlation functional corresponding to 
£tc DA ' ) is constructed from properties of the homoge- 
neous electron gas in the thermodynamic limit (in 
other words, it is the standard LDA functional), the 
functional leading to e s LDA ' ) is based on the homoge- 
neous electron gas confined to the same finite supercell 
as the quantum system under investigation. The latter 
functional is not universal and needs to be found for 
each simulation cell separately at the cost of auxiliary 
simulations of the homogeneous Coulomb gas. 



3.4 An alternative model for Coulomb 
interaction energy 

The expression for the electron-electron interaction 
energy V ee has two properties: (i) it is periodic and 
(ii) corresponds to an actual, albeit artificial, system of 
point charges. Although the latter property is concep- 
tually convenient, it is not really necessary, and any 
periodic potential that exhibits the correct behaviour 
in the limit of the infinitely large simulation cell is le- 
gitimate. Relaxation of the constraint (ii) in favour of 
faster convergence of the total energy per particle 6n 
to its thermodynamic limit was explored in a series 
of papers |113| I120| I125j , where the so-called model 
periodic Coulomb (MPC) interaction was proposed. 
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The replacement for the Ewald energy V t 
1 



(E) 



reads 



v e ( r c) m = 



E 

j<i 



(44) 



E 



V E (ri - r) 



p(r)d 3 



V E (r - r') - 



1 



r — r' 



p(r)p(r')d 3 rd 3 /. 



where \r — r'| m = minR s \r — r'+Rs\ stands for the so- 
called minimum image distance. The operator V r i e MPC ' ) 
is constructed in such a way that the Hartree part of 
its expectation value is the same as in the Ewald 
method, whereas the slowly converging contribution 
to the exchange-correlation energy is removed. There- 
fore, the MPC interaction is essentially equivalent to 
the Ewald formula ( 36 ) complemented with the a pos- 
teriori correction (39). Instead of the structure factor, 
it is the one-particle density p that has to be evaluated 
as an extra quantity in this case (unless it is known 
exactly as in a homogeneous system). The explicit 
presence of the density p in the hamiltonian is incon- 
venient for the DMC method where the local energy 
is needed from the start of the simulation, that is, 
before the density data could be accumulated. The 
situation can be remedied by replacing the unknown 
density p with an approximation p^. For instance, the 
one-particle density provided by DFT is usually quite 
accurate. The error introduced by this substitution is 
proportional to (p — pa) 2 and further diminishes as 
the simulation cell increases. The Ewald and MPC 
energies per particle are therefore the same in the 
thermodynamic limit even if the approximate charge 
density is used PHHH]. 

4 Trial wave functions 

Accurate trial wave functions are essential for success- 
ful applications of the quantum Monte Carlo methods. 
The quality of the employed wave functions influences 
the statistical efficiency of the simulations as well as 
the accuracy of the achieved results. Equally impor- 
tant, especially for investigations of extended systems, 
is the possibility to quickly compute the wave func- 
tion value and its derivatives (V^t and V 2, 3/t), since 
these quantities usually represent the most computa- 
tionally costly part of the whole simulation. Compact 
expressions are therefore strongly preferred. 

A significant part of the construction of the trial 
wave functions is optimization of the variational pa- 



rameters introduced into the functional form repre- 
senting ^t- It is a non-trivial task since the number of 
parameters is often large, <J/t depends non-linearly on 
them, and the quantity to be minimized (£Vmc or the 
variance of the local energy) is a fluctuating number. 
Several powerful methods addressing these problems 
have been developed during the years [127H130] and 
even hundreds of parameters can be optimized with 
good efficiency nowadays. 

4.1 Elementary properties 

Since our aim is the electronic structure, and the elec- 
trons are subject to the Pauli exclusion principle, our 
trial wave functions have to be antisymmetric with re- 
spect to pair-electron exchanges. We assume collinear 
spins that are independent of electron positions, and 
therefore the full wave function can be factorized 
as 



Nl 



c 



x* T (Cft)|C{t^J;l^a}), (45) 

where S = (ci, . . . , <7jv) is a vector consisting of 
N — N-f + N± spin variables. The sum runs over 
all distinct states of N+ up-oriented and Nt down- 



4- 

oriented spins. In the case of = 2 and N± = 1 the 
spin states are | fit2^3>, llYhs^) and |t2T3>k), and 
the corresponding C1Z combinations are {T , 1 ,r 2 ,r 3 }, 
{ri, r 3 , r 2 } and {r 2 , r 3 , ri}. The spatial-only part \f t 
is antisymmetric with respect to exchanges of par- 
allel electrons and its symmetry with respect to ex- 
changes of antiparallel electrons is unrestricted. The 
sum in (45) with the appropriate sign factors (— l) c 
represents the residual antisymmetrization for the an- 
tiparallel spins. 

Both \E't and \&t are normalized to unity and 
identity {9 T \A\^ T ) = (* T |A|* T ) holds for any 
spin-independent operator A since the spin states 
|C{t ■ • • t-l • • • I}) are mutually orthonormal. There- 
fore, it is usually sufficient to consider only the spatial 
part of the full many-body wave function in appli- 
cations of the VMC and DMC methods, and we limit 
our discussion to $t from now on[^] 

Our goal is for the local energy H^t/^t to be 
very close to a hamiltonian eigenvalue and fluctuating 
as little as possible. In systems with charged parti- 
cles interacting via the Coulomb potentials, it requires 
that the kinetic energy proportional to V 2 $t contains 

2 In fact, the DMC algorithm is defined only for the spatial part V^xi consult sections |2. 2. 2| and |2. 2. 4| Decomposition i ]45[ i then 
provides a hint how to calculate expectation values of spin-dependent operators from the sampled mixed distribution \t , J ) \pT' 
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singularities which cancel the 1/r divergencies of the 
potential. This cancellation is vital for controlling the 
statistical uncertainties of the Monte Carlo estimate 
of the total energy and takes place when Kato cusp 
conditions are fulfilled |131| I132j . 

At electron-electron coincidences, these conditions 
can be formulated with the aid of projections of the 
trial wave function v&x onto spherical harmonics Y/ m 
centered at the coincidence point, 

^" m \r tjl r c . m ,,n \ {ri,rj}) 

= J_ f * T (7J)^(^)dny. (46) 

'ij J4tt 

In this definition, the following notation has been in- 
troduced: rij — \rij\ = \ri — rj\ is the electron-electron 
distance, fiy is the spherical angle characterizing ori- 
entation of the vector r^, and r c . m . = (r^ + rj)/2 is 
the position of the center of mass of the electron pair. 
The cusp conditions can then be written as 



lim 



1 



(o.o) dr, 



for unlike spins and 



lim 



(i.' 



(47) 



(48) 



for like spins. The expressions (47) and (48 1 differ 



because is an °dd function with respect to in 
the latter case, which implies = 0. 

Analogous cusps occur in all-electron calculations 
when electrons approach nuclei. Unless the s-wave 
component , 3/; | ? vanishes (for a general discussion 
see |132j ). it can be shown that 



1 <9* 



(0,0) 



(49) 



where ru is the electron-nucleus distance and Zj is 
the charge of the nucleus. 

A convenient functional form that meets the speci- 
fied criteria is a product of an antisymmetric part \I/a 
and a positive symmetric expression exp(— U COTI ) |133] , 



*t(R) = *A(ft) exp[-[/ corr (ft)] 



(50) 



The Jastrow correlation factor exp(- 
porates the electron-electron cusp conditions (47) 



■E/corr) incor- 



and ( 48 ) , and 4»a ensures the fcrmionic character of 



the wave function. The electron-ion cusps ( 49 ) can be 



included in either * A [Ml EU HUH EES] or in the corre- 
lation factor |136j . In simulations of extended systems, 



the antisymmetric part obeys the twisted boundary 
conditions (section 3.1) and exp(— J7 corr ) is periodic 
at the boundaries of the simulation cell. We discuss 



the individual parts of the trial wave function ( 50 ) in 
some detail next. 

4.2 Jastrow factor 

The majority of applications fits into a framework set 
by the expression 



U COII (U) 



j<i 



(erj (r i( rj-), (51) 



where the functions x an d u take a specific 
parametrized form [6lJ I137[ 1138] and can depend also 
on spins of the involved electrons as indicated by the 
indices a e {t>l}- The inclusion of the uncorrelated 
one-body terms \ IS important especially if the trial 
wave function is optimized with a fixed antisymmetric 
part *a [SH I10H H39] . The two-body terms u are 
typically simplified to 



j<i 



3<i,I 



( r ij) 



^■Un/Mr^l), (52) 



where u co is an expression corresponding to a ho- 
mogeneous system and the electron-electron-nucleus 
term w ocn takes into account the differences between 
the two-body correlations in high-density regions near 
nuclei and in low density regions far from them. Spin 
indices were dropped to simplify the notation. The 
tt een contribution can usually be short ranged in the 
\rij\ parameter, whereas the simpler u 00 term is prefer- 
ably long ranged in order to approximate the RPA 
asymptotics (ETl) as closely as allowed by the given 
simulation cell |l0 | 1136] . Of course, limited compu- 
tational resources can (and often do) enforce further 
simplifying compromises. In simulations of homoge- 
neous fermion fluids (electron gas, 3 He), on the other 
hand, even higher order correlations were successfully 
included: three-particle I140H143] as well as four- 
particle |144j . 

4.3 Slater Jastrow wave function 

The simplest antisymmetric form that can be used in 



place of \&a in ( 50 ) is a product of two Slater determi- 
nants, 

xA[ri(r\)...^ N y Ni )} 
= det[^(rj)]det[^(rj)], (53) 
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where and ipfn are single-particle orbitals that cor- 
respond to spin-up respectively spin-down electronic 
states and ipn(rf) in the arguments of the determi- 
nants stands for a N& x N a matrix A G ni . In quantum- 
chemical applications, a common strategy to improve 
upon the Slater wave function is to use a linear com- 
bination of several determinants, 

n- dct W = $> a det[< n (rJ)] det[Vi ro (rj)] . 

a 

(54) 

These multi-determinantal expansions are mostly im- 
practical for simulations of solids, since the number of 
determinants required to describe the wave function to 
some fixed accuracy increases exponentially with the 
system size. One case where multiple determinants are 
vital even in extended systems are fixed-node DMC 
calculations of excitation energies, since adherence 
to proper symmetries is essential for validity of the 
corresponding variational theorem (62j 1145] and the 
trial wave functions displaying the correct symmetry 
are not always representable by a single determinant. 



In these instances, however, the expansions (54) are 
short. 

In simulation cells subject to the twisted boundary 



conditions ( 32 ) specified by a supercell crystal momen- 
tum K , the one-particle orbitals tpm are Bloch waves 
satisfying 



1>Km( r + L <*) =i>Km( r )e 



iK L a 



(55) 



where a = 1,2,3 and m is a band index in the super- 
cell. Since the simulation cell contains several primitive 
cells, the crystal has a higher translational symmetry 



than implied by ( 55 ) and the orbitals can be conve- 
niently re-labeled using m — ¥ (fc,m'), where k and m! 
are a momentum and a band index defined with re- 
spect to the primitive cell. The Bloch waves fulfill 
also 

i>K km >(r + lo t )=r K km>(r)e ik - l '>, (56) 

where l a are lattice vectors defining the primitive cell. 
Assuming that the supercell is built as Let — Tla^a 



with integers n a , the momenta k compatible with ( 55 1 
fall onto a regular mesh 



k = K + 2tt 



J i 



h x U 



m h ■ (h x i 3 ) 

h h*h . 33 h x h 



n 2 h ■ (h x h) n 3 l 3 ■ (h x l 2 ) 



(57) 



with indices j a running from to n a — 1. This set of k 
points corresponds to the Monkhorst-Pack grid [146 
shifted by a vector K. 

A number of strategies have been devised to de- 
termine the optimal one-particle orbitals for use in 



the Slater- J astrow wave functions, which certainly 
differ from the Hartree-Fock orbitals that minimize 
the variational energy only when C/ corr = 0. Ideally, 
the orbitals are parametrized by an expansion in a sat- 
urated basis with the expansion coefficients varied to 
minimize the VMC or DMC total energy. The stochas- 
tic noise and the computational demands of the DMC 
method make the minimization of -Edmc extremely 
inefficient in practice. The VMC optimization of the 
orbitals was successfully performed in atoms and small 
molecules of the first-row atoms |130l 11471 1148] , but 
the method is still too demanding for applications to 
solids. 

To avoid the large number of variational param- 
eters needed to describe the single-particle orbitals, 
another family of methods has been proposed. The 
orbitals in the Slater- Jastrow wave function are found 
as solutions to self-consistent-field equations that rep- 
resent a generalization of the Hartree-Fock theory to 
the presence of the Jastrow correlation factor |139tll49< - 
1151] . These methods were tested in atoms as well as 
in solids within the VMC framework |149] 1152] . Un- 
fortunately, the wave functions derived in this way did 
not lead to lower DMC energies compared to wave 
functions with orbitals from the Hartree-Fock theory 
or from the local density approximation [1491 1153] . It 
is unclear, whether the lack of observed improvements 
in the fermionic nodal surfaces stems from insufficient 
flexibility of the employed correlation factors or from 
the fact that only applications to weakly correlated 
systems were considered so far. 

An even simpler approach is to introduce a para- 
metric dependence into the self-consistent-field equa- 
tions without a direct relation to the actual Jastrow 
factor. An example are the Kohn-Sham equations cor- 
responding to some exchange-correlation functional, in 
which it is possible to identify a parameter (or several 
parameters) measuring the degree of correlations in 
the system and thus mimicking, to a certain degree, 
the effect of the Jastrow factor. Particular hybrid func- 
tional with variable admixture of the exact- exchange 
component |37] were successfully employed for this 
purpose in conjunction with the DMC optimization, 
so that the variations of the fermionic nodal struc- 
ture could be directly quantified [154H157] . Sizeable 
improvements of the DMC total energy associated 
with the replacement of the Hartree-Fock (or LDA) 
orbitals with the orbitals provided by the optimal hy- 
brid functional were observed in compounds containing 
3d elements. 

Evaluation of the Slater determinants dominates 
the computational demands of large-scale Monte Carlo 
calculations, and it is therefore very important to con- 
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sider its implementation carefully. Notably, schemes 
combining a localized basis set (atom-centered Gaus- 
sians or splines |158l [TSQ^ ) with a transformation of 
the single-particle orbitals into localized Wannier-like 
functions can achieve nearly linear scaling of the com- 
putational effort with the system size when applied to 
insulators [T50HT52] . 

4.4 Antisymmetric forms with pair 
correlations 

Apart from the Pauli exclusion principle, the Slater de- 
terminant does not incorporate any correlations among 
the electrons, since it is just an antisymmetrized form 
of a completely factorized function, that is, of a prod- 
uct of one-body orbitals. A better account for cor- 
relations can be achieved by wave functions built as 
the appropriate antisymmetrization of a product of 
two-body orbitals. The resulting antisymmetric forms 
are called pfafhans and can generally be written as 

mom 



in ( 59 ) reduces to the Slater wave function ( 53 ) when 



<J^ faff (ft) = A 



V m,l ' ' m,2 



II Cr 1 " 1 

m— 1 

N-2N P 

x n ^n n « 



(58) 



where Np is the number of correlated pairs, iVp < N/2. 
The two-body orbitals that couple unlike-spin elec- 
trons (singlet pairs) are symmetric, whereas <f)j£ and 
(j>\^ (triplet pairs) are antisymmetric functions. The 
inclusion of the one-body orbitals ip* allows for odd N 
or for only partially paired electrons. The pfaffian 
wave functions can be viewed as compacted forms of 



particular multi-determinantal expansions (54). 



An important representative of the functional 



form (58) is the Bardeen-Cooper-Schrieffer (BCS) 



wave function |164j projected onto a fixed number 
of particles, which is obtained from (58) by con- 



sidering a singlet pairing in an unpolarized system 
(iV-j- = iVj. = N/2) with all two-body orbitals identi- 
cal. In that case, the antisymmetrization reduces to a 
determinant 



(59) 



*l cs (Tl)=det[4>Hrl,rj)], 



where </>^(rJ, rjj) is to be understood as a N/2 x N/2 
matrix Ay. In the quantum-chemical literature, this 
functional form is also known as the antisymmetrized 
geminal power. Note that the form of the BCS wave 
function does not by itself imply formation of Cooper 
pairs and their condensation, since the determinant 



the pair orbital is taken in the form 



y Slator 



N/2 

(rM) = I>U»M(rjo- 

n=l 



(60) 



The BCS-Jastrow wave functions were employed in 



investigations of ultra cold atomic gases (section 5.8) 



1 I167| as well as in calculations of the electronic 
structure of atoms |168j and molecules |169) . 

Trial wave functions with triplet pairing among 
particles were suggested in the context of liquid 3 He 
already two decades ago [1651 1170] . It was realized 
only recently that even in these cases the exponentially 
large number of terms constituting the pfaffian can be 
rearranged in a way that facilitates its evaluation in a 
polynomial time, and therefore allows application of 
the pfaffian-Jastrow trial wave functions in conjunc- 
tion with the VMC and DMC methods [1551 U7I] . 

4.5 Backflow coordinates 

Another way to further increase the variational free- 
dom of the antisymmetric part of the trial wave func- 
tion is the backflow transformation ty^lZ) — > ^a(X), 
where the new collective coordinates X are functions 
of the original electron positions 1Z. The designation 
"backflow" originates from an intuitive picture of the 
correlated motion of particles introduced by Feynman 
to describe excitations in quantum fluids |172l 1173] . 

In order to illustrate what is the origin of 
such coordinates, let us consider homogeneous 
interacting fermions in a periodic box with a 
trial wave function of the Slater- Jastrow type, 
*tM = dct[cxp(ifc • ri)] exp[^ i<3 - l(nj)] ■ The Jas- 
trow factor is optimized so that its laplacian cancels 
out the interactions as much as possible within the 
variational freedom. Applying the kinetic energy op- 
erator to the Slater-Jastrow product results in local 
energy of the form 



[H^ T }(K) 



E var (TZ) 



$ T (K) 

V ln|det [exp(ifc • r*)] 



V^ 7 (r 



i<3 



, (61) 



where we can qualitatively characterize E val {TZ) 
as a mildly varying function close to a con- 
stant while the second term represents a non- 
constant "spurious" contribution, which appears as 
a scalar product of two fluxes. Consider now 
the following modification of the Slater-Jastrow 
form, ^ T (TZ) = det[exp(ife • Xi)] exp[X) i<:j -7(ry)], 
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where the single-particle coordinates are modified as 
Xi = fi+J^j ^( r ij) with #(0) = 0. One can show that 
with a proper choice of the function i9(r^), the lapla- 
cian of det [exp(ifc • asj)] produces terms that cancel out 
most of the spurious contributions in the local energy 



given by (61). Of course, the backflow form generates 



also new non-constant terms so that the wave function 
has to be optimized for the overall maximum gain 
using variational strategies. 

In general, the new coordinates are written as 
Xi = ri + £j (1Z) with £ taken in a form analogous to 
the parametrization of the Jastrow factor U con - (51) 



and ( 52 ) . The vector £ contains two-particle and pos- 



sibly higher order correlations and, in systems with ex- 
ternal potentials, also inhomogeneous one-body terms. 
The backflow transformation has been very success- 
ful in simulations and understanding of homogeneous 
quantum liquids 11411 11431 1144j , and some progress 
has recently been reported in applying these techniques 
also to atoms and molecules |163l 1174) . 

5 Applications 

In the last part of the article we go through selected 
applications of the quantum Monte Carlo methodology 
to the electronic structure of solids. In practically all 
listed cases, with the exception of sections [5~T] and [578] 
dealing with model calculations, the Slater- Jastrow 
functional form is employed as the trial wave function. 
The reviewed results therefore map out the accuracy 
that is achievable in the realistic solids when the mean- 
field topology of the fermionic nodes is assumed. It is 
shown that the quality of the DMC predictions is re- 
markable despite this relatively simple approximation. 

5.1 Properties of the homogeneous electron 
gas 

The homogeneous electron gas, also referred to as jel- 
lium, is one of the simplest many-body models that 
can describe certain properties of real solids, especially 
the alkali metals. At zero temperature, the model 
is characterized by the densities of spin-up and spin- 
down electrons, p^ and pi, or, alternatively, by the 
total density p = p^ + p± and the spin polarization 
C = — Pi\/p- It is convenient to express the den- 
sity p and other quantities in terms of a dimensionless 
parameter r s = [3/(47r y o)] 1/ ' 3 /aB, where ae is the Bohr 
radius. For example, the density of valence electrons 
in the sodium metal corresponds to r s « 4. 

The total energy of jellium is particularly simple 
since it includes only the kinetic energy of the electrons, 
the Coulomb electron-electron repulsion, and a con- 
stant which represents the interaction of the electrons 



with an inert uniformly distributed positive charge 
that maintains overall charge neutrality of the system. 
A straightforward dimensional analysis shows that the 
kinetic energy dominates the Coulomb interaction at 
high densities (small r s ), where the electrons behave 
like a nearly ideal gas and the unpolarized state (C = 0) 
is the most stable. In the limit of very low densities, on 
the other hand, the kinetic energy becomes negligible 
and the electrons "freeze" into a Wigner crystal [189J. 

The homogeneous electron gas at zero temperature 
was one of the first applications of the variational and 
diffusion Monte Carlo methods. In the early inves- 
tigations [ini [H], only the unpolarized (C = 0) and 
fully polarized (C = 1) fluid phases were considered 
together with the Wigner crystal. Later, fluids with 
partial spin polarization were taken into account as 
well [190 193 . The most accurate trial wave functions 
(the Slater-Jastrow form with backflow correlations) 
were used in reference |193j where it was found that 
the unpolarized fluid is stable below r s = 50 ± 2. At 
this density the gas undergoes a second-order phase 
transition into a partially polarized state, and the 
spin polarization £ then monotonically increases as 
the fluid is further diluted. Eventually, the Wigner 
crystallization density is reached, for which two DMC 
estimates exist: r s = 100 ± 20 [llj and r s — 65 ± 10 



[192 . The discrepancy is presumably caused by the 
very small energy differences between the competing 
phases over a wide range of densities, and by uncer- 
tainties in the extrapolation to the thermodynamic 
limit. Advanced finite-size extrapolation methods, out- 
lined in section [3] earlier, could possibly shed some new 
light on these quantitative differences. Indeed, recent 
calculations show further improvements in accuracy 
of the total and correlation energies |194j . A number 
of static properties of the liquid phases that provide a 
valuable insight into the details of the electron corre- 
lations in the jellium model and in Coulomb systems 
in general were evaluated by QMC methods as well 

The impact of the QMC calculations of the ho- 
mogeneous electron gas |llj has been very significant 
because of the prominent position of the model as one 
of the simplest periodic many-body systems, and also 
due to the fact that the QMC correlation energy has 
become widely used as an input in density-functional 
calculations [3"5l 1196j . 

The results quoted so far referred to the homo- 
geneous Coulomb gas in three dimensions. The two- 
dimensional gas, which is realized by confining elec- 
trons to a surface, interface or to a thin layer in a 
semiconductor heterostructure, has received similar 
if not even greater attention of QMC practitioners 
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Table 1: Cohesive energies of solids (in eV). Shown are 
DMC numbers unless only VMC data are available for 
the particular compound; those instances are marked 
with (*). The latest results are preferred in cases where 
multiple calculations exist. If not indicated otherwise, 
the experimental cohesive energies are deduced from the 
room-temperature formation enthalpies quoted in |188j . 



compound 


QMC 




experiment 


t : 

J_j1 


1 HQ -1- n n^i 

i.uy zt u.uo 
1.57±0.01 


1 ' o. 

[T76] (*) 


1.00 


Na 


1.14±0.01 
1.0221 ±0.0003 


[124] 
[126] 


l.n 


Mg 


1.51±0.01 


EH 


1.52 


Al 


3.23 ±0.08 


HZ3 (*) 


3.43 [HZ] 


MgH 2 


6.84±0.01 


una 


6.83 


BN 


12.85 ±0.09 


mi (*) 


12.9 [T79] 


C (diamond) 


7.346 ±0.006 


mi 


7.37 [HI] 


Si 


4.62 ±0.01 


[152] 


4.62 [THS] 


Gc 


3.85 ±0.10 


mi 


3.86 


GaAs 


4.9 ±0.2 


nn (*) 


6.7 [US] 


MnO 


9.29 ±0.04 


EST] 


9.5 


FcO 


9.66 ±0.04 


m 


9.7 


NiO 


9.442 ±0.002 


mi 


9.5 


BaTiOs 


31.2±0.3 


nsa 


31.57 



[TU1 [T?7rf2U2] . The Wigner crystallization was pre- 
dicted to occur at r s = 37 ± 5 l~l!>7 a value that 
coincides with the density r s — 35 ± 1 where a metal- 
insulator transition was experimentally observed later 
[203] , 

5.2 Cohesive energies of solids 



The cohesive energy measures the strength of the chem- 
ical bonds holding the crystal together. It is defined as 
the difference between the energy of a dilute gas of the 
constituent atoms or molecules and the energy of the 
solid. Calculation of the cohesive energy is a stringent 
test of the theory, since it has to accurately describe 
two different systems with very dissimilar electronic 
structure. 

The first real solids whose cohesive energies were 
evaluated by a QMC method were carbon and silicon 
in the diamond crystal structure [1011 1204] . These 
early VMC estimates were later refined with the DMC 
method [MB EEE1 HE3 HOI]. The most accurate re- 
sults to date are shown in table [I] where we have 
compiled the cohesive energies of a number of com- 
pounds investigated with the quantum Monte Carlo 
methods. Corresponding experimental data are shown 
for comparison. The electronic total energy calculated 
in QMC simulations is not the only contribution to 
the cohesive energy of a crystal, and the zero-point 
and thermal motion of the nuclei has to be accounted 
for as well, especially in compounds containing light 

3 Note that in two dimensions the dimensionless density parameter r a is defined as r. 



atoms. We refer the reader to the original references 
for details of these corrections. At present, a direct 
QMC determination of the phonon spectrum is gen- 
erally not practicable due to unresolved issues with 
reliable and efficient calculation of forces acting on the 
nuclei [IB]. The effects due to the nuclear motion are 
thus typically estimated within the density-functional 
theory. 



Overall, the agreement of the DMC results with 
experiments is excellent; the errors are smaller than 
0.1 eV most of the time, including the Na and Mg 
elemental metals where coping with the finite-size ef- 
fects is more difficult. Notably, the diffusion Monte 
Carlo performs (almost) equally well in strongly cor- 
related solids represented in table [l] by 3d transition 
metal oxides MnO, FeO and NiO. The GaAs result is 
an obvious outlier with a systematic error of almost 
2 eV that the authors identify with the deficiencies 
of their pseudopotentials [184J. The two decades old 
application of the DMC method to the Li metal [175] 
is the only all-electron simulation in the list and its 
comparison to a subsequent pseudopotential calcula- 
tion [176] suggests that a large part of the discrepancy 
with the experiment is due to the fixed-node errors 
in the high-density core regions. It is likely that a 
substantial improvement would be observed if the all- 
electron calculations were revisited with today's state 
of the art trial wave functions. 

l/(a Bv / 7fp)- 
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compound 


ao (A) 




Vb (A 3 ) 




Bo (GPa) 




Li 


3.556 ± 0.005 
3.477 


[176] 
[207] 






13 ±2 
12.8 


[126] 
[205] 


Al 


3.970 ±0.014 [177] 
4.022 [T77] 






65 ± 17 
81.3 


[177] 
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GaAs 


5.66 ±0.05 
5.642 


.184 
[209] 
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77 ±1 


ma 

[2T0] 


LiH 


4.006 

4.07 ±0.01 


[2TT] 

ma 
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32.2 ±0.03 


[2TT] 
EH] 


BN 






11.812 ±0.008 
11.812 ±0.001 


[T55] 
[2T3] 


378 ±3 
395 ±2 
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[2T3] 


Mg 






23.61 ±0.04 
23.24 


106 
214 


31.2 ±2.4 
36.8 ±3.0 


[TUB] 
215 


MgO 


4.23 
4.213 


ms] 
[nz] 






158 

160 ± 2 


[216] 

217 


MgH 2 






30.58 ±0.06 
30.49 


[106] 

EH 


39.5 ± 1.7 


ESI 


C 

(diamond) 


3.575 ±0.002 
3.567 


[2T9] 

BBS1 






437 ±3 
442 ±4 


[2T5] 
[220] 


Si 


5.439 ± 0.003 [152] 
5.430 [221] 






103 ± 10 
99.2 


HFJ 
[222] 


SiOa 
(quartz) 






37.6 ±0.3 
37.69 


[223] 
[224] 


32 ±6 
34 


[223] 
[224] 


FeO 


4.324 ± 0.006 [TTS] 
4.334 [225] 






170 ± 10 
w 180 


.119 
[226] 



Table 2: Equilibrium lattice constants 
Oo, equilibrium volumes Vo (per formula 
unit) and bulk moduli Bo for a number of 
solids investigated with the QMC meth- 
ods. The first line for each compound 
contains QMC predictions, the second 
line shows experimental data. Theoret- 
ical results for Li, Al and GaAs come 
from VMC simulations, the rest of the 
table corresponds to the DMC method. 



5.3 Equations of state 

The equilibrium volume Vo, the lattice constant do, and 
the bulk modulus B = V(dE /dV)\v=v are among 
the most basic parameters characterizing elastic prop- 
erties of a solid near the ambient conditions. Within 
QMC methods, these quantities are determined by eval- 
uating the total energy at several volumes around Vq 
and by fitting an appropriate model [2271 1228j of the 
equation of state E(V) through the acquired data, 
see figure [2] for an illustration. Results of this proce- 
dure for a wide range of solids are shown in table [2] 
together with the corresponding experimental data. 
As in the case of the cohesion energy discussed in the 
preceding section, the raw QMC numbers correspond 
to the static lattice and corrections due to the motion 
of nuclei may be needed to facilitate the comparison 
with experimentally measured quantities. Particular 
details about applied adjustments can be inspected in 
the original papers. 

The data in table [2] demonstrate that the equilib- 
rium geometries predicted by the QMC simulations are 



very good and all lie within 2 % from the experiments, 
in many cases within only a few tenths of a percent. 
The agreement is slightly worse for the bulk moduli, 
where errors of several percent are common and in 
a few instances the mean values of the Monte Carlo 
estimates deviate from the experimental numbers by 
more than 10 %. Note, however, that determination of 
the curvature of E(V) near its minimum is impeded by 
the stochastic noise of the QMC energies and that the 
error bars on the less favourable results are relatively 
large. 

Quantum Monte Carlo methods are not limited 
to the covalent solids listed in table [2| Investigation 
of the equation of state of solid neon [229] represents 
an application to a crystal bound by van der Waals 
forces. Although the shallowness of the minimum 
of the energy-volume curve in combination with the 
Monte Carlo noise did not allow to determine the 
lattice constant and the bulk modulus to a sufficient 
accuracy, the DMC equation of state was still substan- 
tially better than results obtained with LDA and GGA. 
This example together with a recent study of interlayer 
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Figure 2: DMC total energies of the rock-salt 
(squares) and the NiAs (circles) phases of FcO. 
Statistical error bars are smaller than the symbol 
sizes. Lines are fits with the Murnaghan equa- 
tion of state. Inset: Difference between the Gibbs 
potentials of the two phases at T = K; its in- 
tercept with the x axis determines the transition 
pressure Pf Adopted from |119| . 




14 16 18 20 22 24 

V(A 3 /FeO) 



binding in graphite |230j illustrate that the diffusion 
Monte Carlo method provides a fair description of dis- 
persive forces already with the simple nodal structure 
defined by the single-determinantal Slater- Jastrow 
wave function. More accurate trial wave functions 
incorporating backffow correlations were employed to 
study van der Waals interactions between idealized 
metallic sheets and wires [231] . 

Calculations of the equations of state are by no 
means restricted to the vicinity of the equilibrium vol- 
ume, and many of the references quoted in table [2] 
study the materials up to very high pressures. Such 
investigations are stimulated by open problems from 
Earth and planetary science as well as from other areas 
of materials physics. Combination of the equation of 
state with the pressure dependence of the Raman fre- 
quency |135l 1219] , both calculated from first principles 
with QMC, can provide a very accurate high-pressure 
calibration scale for use in experimental studies of 
condensed matter under extreme conditions [135J. 

5.4 Phase transitions 

Theoretical understanding of structural phase transi- 
tions often necessitates a highly accurate description 
of the involved crystalline phases. Simple approxi- 
mations are known to markedly fail in a number of 
instances due to significant changes in the bonding 
conditions across the transition. A classical example is 
the quartz-stishovite transition in silica (Si02), where 
LDA performs very poorly and GGA is needed to 
reconcile the DFT picture with experimental findings 
[232J . The diffusion Monte Carlo method has been 
employed to investigate pressure-induced phase tran- 



sitions in Si PH], MgO [US], FeO [112] and Si0 2 
[22~3] . 

A transition from the diamond structure to the 
/3-tin phase in silicon was estimated to occur at 
P% = 16.5 ±0.5 GPa [182] . which lies outside the range 
of experimentally determined values 10.3-12.5 GPa 
(see 182 j for compilation of experimental literature). 
Since the diamond structure is described very accu- 
rately with the DMC method as testified by the data 
in tables [T] and [2j it was suggested that the discrep- 
ancy is a manifestation of the fixed-node errors in the 
high-pressure /3-tin phase. This view is supported by 
a recent calculation utilizing the so-called phaseless 
auxiliary-field QMC, a projector Monte Carlo method 
that shows smaller biases related to the fermion sign 
problem in this particular case and predicts the tran- 
sition at 12.6 ± 0.3 GPa [233 . It should be noted, 
however, that the volume at which the transition oc- 
curs was fixed to its experimental value in this later 
study, whereas the approach pursued in [182] was en- 
tirely parameter-free. 

In iron oxide (FeO), a transition from the rock- 
salt structure to a NiAs-based phase is experimentally 
observed to occur around 70 GPa at elevated tempera- 
tures |234j and to move to higher pressures exceeding 
100 GPa when the temperature is lowered |235j . The 
DMC simulations summarized in figure [2] place the 
transition at P t = 65 ± 5 GPa [119] . This value repre- 
sents a significant improvement over LDA and GGA 
that both find the NiAs structure more stable than 
the rock-salt phase at all volumes. The agreement 
with experiments is nevertheless not entirely satisfac- 
tory, since the DMC prediction corresponds to zero 
temperature where experimental observations suggest 
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stabilization of the rock-salt structure to higher pres- 
sures. Sizeable sensitivity of the transition pressure Pt 
to the choice of the one-particle orbitals in the Slater- 
Jastrow trial wave function was demonstrated in a 
subsequent study |157] . but those wave functions that 
provided higher P t also increased the total energies, 
and therefore represented poorer approximations of the 
electronic ground state. It remains to be determined, 
whether the discrepancy between the experiments and 
the DMC simulations is due to inaccuracies of the 
Slater- Jastrow nodes or if some physics not included 
in the investigation, for instance the inherently defec- 
tive nature of the real FeO crystals, plays a significant 
role. 

Investigations of phase transitions involving a liq- 
uid phase, such as melting, are considerably more 
involved due to a non-trivial motion of ions. An often 
pursued approach is a molecular dynamics simulation 
of ions subject to forces derived from the electronic 
ground state that is usually approximated within the 
density-functional theory. More accurate results would 
be achieved if the forces were calculated using quan- 
tum Monte Carlo methods instead. At present, this 
is generally not feasible due to excessive noise of the 
available force estimates |48j . Nevertheless, it was 
demonstrated that one can obtain an improved pic- 
ture of the energetics of the simulated system when 
its electronic energy is evaluated with the aid of a 
QMC method while still following the ion trajectories 
provided by the DFT forces [23611237] , 

5.5 Lattice defects 

The energetics of point defects substantially influences 
high-temperature properties of crystalline materials. 
Experimental investigations of the involved processes 
are relatively difficult, and it would be very helpful if 
the electronic structure theory could provide depend- 
able predictions. The role of electron correlations in 
point defects was investigated with the DMC method 
in silicon 206, 238] and in diamond |180j . The forma- 
tion energies of selected self-interstitials in silicon were 
found about 1.5 eV larger than in LDA, whereas the 
formation energy of vacancies in diamond came out 
as approximately 1 eV smaller than in LDA. These 
differences represent 20-30 % of the formation ener- 
gies and indicate that an improved account of electron 
correlations is necessary for accurate quantitative un- 
derstanding of these phenomena. 

Charged vacancies constituting the Schottky defect 
were investigated in MgO [239], and in this case the 
predictions of the DMC method differ only marginally 
from the results obtained within the local density 
approximation. The non-zero net charge of the su- 



percells employed in these simulations represents an 
additional technical challenge in the form of increased 
finite-size effects that require a modification of some of 
the size extrapolation techniques discussed in section [3] 
[240 1 124T ] . 

5.6 Surface phenomena 

Materials surfaces are fascinating systems from the 
point of view of electronic structure and correlation 
effects. The vacuum boundary condition provides 
surface atoms with more space to relax their posi- 
tions and surface electronic states enable the electronic 
structure to develop features which cannot form in 
the periodic bulk. This leads to a plethora of sur- 
face reconstruction possibilities with perhaps the most 
studied paradigmatic case of 7 x 7 Si(lll) surface 
reconstruction. Seemingly, QMC methods should be 
straightforward to apply to these systems, similarly to 
the three-dimensional periodic solids. However, mainly 
technical reasons make such calculations quite difficult. 
There are basically two possibilities how to model a 
surface. One option is to use a two-dimensional pe- 
riodic slab geometry which requires certain minimal 
slab thickness in order to accurately represent the bulk 
environment for the surface layers on both sides. The 
resulting simulation cells end up quite large making 
many such simulations out of reach at present. The 
other option is to use a cluster with appropriate ter- 
mination that mimics the bonding pattern of the bulk 
atoms. This strategy assumes that the termination 
does not affect the surface geometries in a substantial 
manner. Moreover, it is applicable only to insulating 
systems. Given these difficulties, the QMC simulations 
of surfaces are rare and this research area awaits to 
be explored in future. 

The simplest possible model for investigation of 
surface physics is the surface of the homogeneous elec- 
tron gas that has been studied by DFT as well as QMC 
methods. The first QMC calculations [243] were later 
found to be biased due to complications arising from 
finite-size effects, especially due to different scaling of 
finite-size corrections for bulk and surface. Once these 
issues have been properly taken into account by Wood 
and coworkers |242) . the QMC results have exhibited 
trends that were consistent with DFT and RPA meth- 
ods which are expected to perform reasonably well for 
this model system (see table [3| . 

Applications to real materials surfaces are still very 
few. The cluster model was used in calculations of 
Si(001) surface by Healy et al. |244| with the goal 
of elucidating a long-standing puzzle in reconstruc- 
tion geometry of this surface, which exhibits regularly 
spaced rows of Si-Si dimers. The dimers could take 
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Table 3: Comparison of the surface energies (in erg cm -2 ) of the 
homogeneous electron gas calculated by a number of electronic struc- 
ture methods [242]. The DMC calculations were done with the LDA 
orbitals in the trial wave functions. 



r a 


LDA 


GGA 


DMC 


RPA 


2.07 


-608.2 


-690.6 


-563 ± 45 


-517 


2.30 


-104.0 


-164.1 


-82 ±27 


-34 


2.66 


170.6 


133.0 


179 ± 13 


216 


3.25 


221.0 


201.2 


216 ±8 


248 


3.94 


168.4 


158.1 


175 ±8 


182 



two possible conformations: They can be either posi- 
tioned symmetrically or form an alternating buckling 
in a zig-zag fashion. While experiments suggested the 
buckled geometry as the low-temperature ground state, 
theoretical calculations produced conflicting results, 
in which various methods favoured one or the other 
structure. The QMC calculations |244j concluded that 
the buckled geometry is lower by about 0.2 eV/Si2- 
This problem was studied with QMC methodology 
also by Bokes and coworkers [245] who found that 
several systematic errors (such as uncertainty of ge- 
ometries in cluster models and pseudopotential biases) 
added to about 0.2 eV, and therefore prevented un- 
equivocal determination of the most stable geometry. 
This conclusion corroborated the experimental find- 
ings which suggested that at temperatures above 100 
K the distinct features of buckling were largely washed 
out and indicated that the effect is energetically very 
small. Very recently, the QMC study of this system 
has been repeated by Jordan and coworkers [246] with 
the conclusion that the buckled structure is lower by 
about 0.1 eV/Si2 and that the highest level correlated 
basis set method which they used (CASPT3) is consis- 
tent with this finding. It was also clear that once the 
correlations were taken into account, the energy differ- 
ences between the competing surface reconstruction 
patterns were becoming very small. This brings the 
calculations closer to reality, where the two structures 
could be within a fraction of 0.1 eV/Si2 as suggested 
by experiments. 

Perhaps the most realistic QMC calculations of 
surfaces have been done on LiH and MgO surfaces by 
the group of Alfe and Gillan [2111 1247j who compared 
predictions of several DFT functionals with the fixed- 
node DMC method. The results showed significant 
differences between various exchange-correlation func- 
tionals. For the MgO(100) surface the best agreement 
with QMC results was found for the LDA functional, 
while for the LiH surface the closest agreement between 
QMC and DFT predictions was found for particular 
GGA functionals. 

Clearly, more applications are needed to assess the 
effectiveness of QMC approaches for investigation of 
surface physics. As we have already mentioned, the 



surfaces represent quite challenging systems for QMC 
methods. Nevertheless, we expect more applications 
to appear in the future since the field is very rich in va- 
riety of correlation effects that are difficult to capture 
by more conventional methods. 

5.7 Excited states 

The VMC and the fixed-node DMC methods both 
build on the variational principle, and they therefore 
seem to be applicable exclusively to the ground-state 
properties. Nevertheless, the variational principle can 
be symmetry constrained, in which case the algorithms 
search for the lowest lying eigenstate within the given 
symmetry class (provided, in the case of the DMC 
method, that the eigenstate is non-degenerate [62] h 
and thus enable access to selected excited states. 

Excitation energies in solids are calculated as differ- 
ences between the total energy obtained for the ground 
state and for the excited state. It is a computationally 
demanding procedure since the stochastic fluctuations 
of the total energies are proportional to the number of 
electrons in the simulation cell, whereas the excitation 
energy is an intensive quantity. Trial wave functions 
for excited states are formed by modifying the deter- 
minantal part of the ground-state Slater- Jastrow wave 
function such that an occupied orbital in the ground- 
state determinant is replaced by a virtual orbital. This 
substitution corresponds to an optical absorption ex- 
periment where an electron is excited from the valence 
band into the conduction band. The fact that both the 
original occupied orbital and the new virtual orbital 
necessarily belong to the same K point restricts the 
types of excitations that can be studied, since only a 
limited number of k points from the primitive cell fold 
to the given K point of the simulation cell, recall equa- 



tions (55|-(57|. Clearly, the larger the simulation cell, 



the finer mapping out of excitations can be performed. 

Averaging over twisted boundary conditions (sec- 
tion 3.1) is not applicable to the calculations of the 



excitation energies, since both the ground state and 
the excited state are fixed to a single K point. This 
is not a significant issue, since finite-size effects tend 
to cancel very efficiently in the differences of the total 
energies calculated at the same K point. 
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Table 4: Band gaps (in eV) of Mott insulators MnO and FeO cal- 

compound PMC experiment culated with the fixed _ node DMC met hod. Experimental data are 

MnO 4.8 ± 0.2 [248] 3.9 ± 0.4 [249] provided for comparison. 

FeO 2.8 ± 0.3 [119] « 2.4 [250] 



DMC simulations following the outlined recipe were 
utilized to estimate the band gap in solid nitrogen 
[25 Ij and in transition-metal oxides FeO [119] and 
MnO |248j . The gaps calculated for the two strongly 
correlated oxides are compared with experimental data 
in table [4] The ratio of the FeO and MnO gaps is 
reproduced quite well, but the DMC gaps themselves 
are somewhat overestimated, likely due to inaccuracies 
of the trial wave functions used for the excited states. 
A large number of excitations were calculated in sili- 
con |252j and in diamond |1451 [253] , and the obtained 
data were used to map, albeit sparsely, the entire 
band structure. In these weakly correlated solids the 
agreement with experiments is very good. Recently, a 
pressure-induced insulator-metal transition was inves- 
tigated in solid helium by calculating the evolution of 
the band gap with compression [254 . As illustrated in 
figure [3] (copyrighted material; not available in this ver- 
sion; see figure 1 in \25^j ), the DMC band gaps were 
found to practically coincide with the gaps calculated 
with the GW method. 



A large amount of research activity aimed at detailed 
understanding of this physics was stimulated by the 
possibility to realize the BCS-BEC crossover in exper- 
iments with optically trapped ultra-cold atoms [259] . 

In a dilute Fermi gas with short-ranged spherically 
symmetric inter-particle potentials, the interactions 
are fully characterized by a single parameter, the two- 
body scattering length a. The system is interpolated 
from the BCS regime to the BEC limit by varying 
1/a from — 00 to 00. In experiments, this is achieved 
by tuning across the Feshbach resonance with the aid 
of an external magnetic field. Particularly intriguing 
is the quantum state of an unpolarized homogeneous 
gas at the resonance itself, where the scattering length 
diverges (1/a = 0). The only relevant length scale 
remaining in the problem in this case is the inverse 
of the Fermi wave vector 1/fcp, and all ground-state 
properties should therefore be universal functions of 
the Fermi energy E-p. Since there is only a single 
length scale, the system is said to be in the unitary 
limit. The total energy can be written as 



5.8 BCS BEC crossover 

The repulsive Coulomb interaction considered so far 
is not the only source of non-trivial many-body ef- 
fects in the electronic structure of solids. A weak 
attractive interaction among electrons is responsible 
for a very fundamental phenomenon — the electronic 
states in the vicinity of the Fermi level rearrange into 
bosonic Cooper pairs that condense and give rise to 
superconductivity. The ground state of the system 
can be described by the BCS wave function ^bcs dis- 

The Cooper pair 



E = SEf, 



cussed earlier in section 4.4 



is an entity that has a meaning only as a constituent 
of "Sbcs- In order to form an isolated two-electron 
bound state, some minimal strength of the two-body 
potential is needed in three dimensions, whereas the 
Cooper instability itself occurs for arbitrarily weak 
attraction. When the interaction is very strong, the 
composite bosons formed as the two-electron bound 
states indeed exist and undergo the Bose-Einstein 
condensation (BEC). It turns out that a mean-field 
description of both the BCS and BEC limits leads to 
the same form of the many-body wave function, which 
indicates that the interacting fermionic system is likely 
to continuously evolve from one limit to the other when 
the interaction strength is gradually changed 256 258 . 
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(62) 



where Ef rcc denotes the energy of a non-interacting 
system and £ is a universal parameter. The universal- 
ity of £ is illustrated in figure [4] that shows the ratio 
E I Bfree as a function of the interaction strength cal- 
culated for three different particle densities using the 
diffusion Monte Carlo method with the trial wave func- 
tion of the BCS-Jastrow form. All three curves indeed 
intersect at 1/a = with the parameter £ estimated 
as 0.42 ± 0.01 [TB5II2WJ352] . The energy calculated 
with the fermionic nodes fixed by the Slater- Jastrow 
wave function is considerably higher and would lead 
to £ ps 0.54 [166] . which underlines the significance of 
particle pairing in this system. 

A further insight into the formation of the Cooper 
pairs is provided by evaluation of the condensate frac- 
tion that can be estimated from the off-diagonal long- 
range order occurring in the two-particle density ma- 
trix |263j . The condensate fraction a is given as 



N 



2 r-y 



lim Pz n (r) 



(63) 



and the so-called projected two-particle density matrix 
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Figure 4: Fixed-node DMC energies of 38 
fermions in a cubic box with the periodic bound- 
ary conditions plotted as a function of interaction 
strength. BEC regime is on the left, BCS limit 
on the right. Shown are three particle densities p 
characterized by the dimensionless parameter r s 
defined in section [57T] The simulations employed 
BCS-Jastrow trial wave function. Statistical er- 
ror bars are smaller than the symbol sizes. Data 
taken from 255 . 
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where r± corresponds to the spin-up state and r 2 to 
the spin-down state. The evolution of a with interac- 
tion strength calculated with the DMC method [265 
is shown in figure [5] ( copyrighted material; not avail- 
able in this version; see figure 4 in 12651). It is found 
that approximately half of the particles participates in 
pairing in the unitary regime and this fraction quickly 
decreases towards the BCS limit, where only states in 
the immediate vicinity of the Fermi level contribute to 
the Cooper pair formation. Note that the condensate 
fraction vanishes if the Slater-Jastrow form is used in 
place of the trial wave function. 

The diffusion Monte Carlo simulations were used 
to study also the total energy and the particle density 
profile in the unitary Fermi gas subject to a harmonic 
confining potential |266H268j . Due to the lowered 
symmetry compared to the homogeneous calculations 
referred above, the system sizes were more limited. To 
extrapolate the findings to a larger number of particles, 
a density functional theory fitted to the DMC data 
can be employed |269j . 



6 Concluding remarks 

In this article we have attempted to provide an 
overview of selected quantum Monte Carlo methods 
that facilitate calculation of various properties of cor- 
related quantum systems to a very high accuracy. Par- 
ticular attention has been paid to technical details 
pertaining to applications of the methodology to ex- 
tended systems such as bulk solids. We hope that we 
have been able to demonstrate that the QMC methods, 
thanks to their accuracy and a wide range of applica- 
bility, represent a powerful and valuable alternative to 
more traditional ab initio computational tools. 
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